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Abstract 

In a recent paper pQ we presented precise lattice QCD results of our Eu- 
ropean Twisted Mass Collaboration (ETMC). They were obtained by em- 
ploying two mass-degenerate flavours of twisted mass fermions at maximal 
twist. In the present paper we give details on our simulations and the com- 
putation of physical observables. In particular, we discuss the problem of 
tuning to maximal twist, the techniques we have used to compute corre- 
lators and error estimates. In addition, we provide more information on 
the algorithm used, the autocorrelation times and scale determination, the 
evaluation of disconnected contributions and the description of our data 
by means of chiral perturbation theory formulae. 
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1 Twisted mass fermions 

Dynamical Wilson twisted mass fermions, when tuned to maximal twist [21 E], 
have been demonstrated to lead to precise results for mesonic quantities down 
to pseudo scalar masses mps < 300 MeV. Results in the quenched case were 
discussed in refs. [H El E] and in the case of two mass-degenerate flavours of 
quarks in ref. pQ. Preparatory simulations with twisted mass dynamical fermions 
were performed in [TJElElEin]- hi ref. pQ many of the details of our computations 
had to be omitted and it is the purpose of the present paper to supplement those 
and fill this gap. 

This paper is organized as follows. In this section we introduce twisted mass 
fermions and discuss the important issue of tuning to maximal twist. In section [21 
we give details about our techniques to compute charged correlators and in sec- 
tion El to compute neutral correlators and quark-disconnected contributions. In 
section H] we discuss the algorithm details and explain our analysis techniques for 
obtaining reliable error estimates. In section E] we provide details of our com- 
putation of the force parameter tq and in section [6] we give some results for the 
pseudoscalar mass and decay constant, the untwisted PCAC quark mass and the 
renormalization constant Zy. We use chiral perturbation theory to fit our data 
and we detail this procedure in section [71 We end with a short summary in 
section El 

We begin with the Wilson twisted mass fermionic lattice action for two 
flavours of mass degenerate quarks, which reads (in the so called twisted ba- 
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sis [21 E] and with fermion fields having continuum dimensions) 

= fl4 E{^ [Dvf + m + Z7 5 r 3 /i 9 ] x*} , 

(1) 

d w = -7, (v, + v; t ) - - v,v; , 

where mo is the bare untwisted quark mass and \x q the bare twisted quark mass, 
r 3 is the third Pauli matrix acting in flavour space and r is the Wilson parameter, 
which we set to r = 1 in our simulations. We denote by and V* the gauge 
covariant nearest neighbour forward and backward lattice derivatives. The bare 
quark mass mo is related to the so-called hopping parameter k, which we will 
often use in this paper, by k = 1/(8 + 2am ). Twisted mass fermions are said 
to be at maximal twist if the bare untwisted mass is tuned to its critical value, 
m crit- We will discuss later how this can be achieved in practice. 

In the gauge sector we use, for reasons explained in pQ, the so called tree- 
level Symanzik improved gauge action (tlSym) [II] which includes besides the 
plaquette term 17** \, also rectangular (1x2) Wilson loops 17**^, ■ It reads 

s * = § E ( 6 ° E i 1 - ReTr (^)}+ & i E o - R eTr (C)} ) > ( 2 ) 

i<fj,<v i±¥ zv 

where (3 is the bare inverse coupling and we set b\ = —1/12 (with b = 1 — 8&i 
as dictated by the requirement of continuum limit normalization). Note that at 
bi = this action becomes the usual Wilson plaquette gauge action. 

1.1 Tuning to maximal twist 

One of the main virtues of Wilson twisted mass fermions is that by tuning the 
bare quark mass mo to its critical value an automatic 0(a) improvement can 
be achieved such that expectation values of parity even operators scale to their 
continuum limit with 0(a 2 ) discretization errors [3]. It was shown in the scaling 
test study carried out in [H El [H] in the quenched case that 0(a) improvement 
works extremely well for maximally twisted mass quarks. In this context, the 
method to tune to maximal twist by setting the so-called (untwisted) PCAC 
mass to zero (in the limit \x q — > 0) was found to be very successful, in agreement 
with theoretical considerations [121 [El E] • In the present paper essentially the 
same approach to set to zero the (untwisted) PCAC mass 

m _ Ex(^g(x,t)P a (o)) 

was followed, by evaluating (jHJ) at large enough time separation, so that the 
pion ground state is dominant. For a definition of the (twisted basis) operators 
appearing in eq. © see eq. (H7I) of Appendix [A] 
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In principle one could think of determining am cr \ t at each value of a/x 9 at 
which simulations are carried out and then perform an extrapolation to vanish- 
ing afi q based on data satisfying the bound a\i q > a 3 AQ CD [H]. This method 
is, however, rather CPU-time expensive. We therefore prefer to determine the 
value of am cr it (at each fixed value of 0) from the simulation at the lowest avail- 
able value, a/ig imin <C aAq CD . This choice simply affects the critical quark mass 
by O (a/i 9 min Aq CD ) terms. Therefore 0(a) improvement is still guaranteed [3J. 
Furthermore, and most importantly, with such a determination of am crit also the 
0(a 2 ) cutoff effects remain small as long as \i q > cl 2 Aq CD [TJ]. We recall below 
the line of arguments leading to this conclusion. 

1.2 Maximal twist and residual 0(a 2 ) artifacts 

To start the discussion let us assume that mo = l/(2n) — 4 has been set to a 
value corresponding to some sensible lattice estimate of the critical mass, while 
/iq is non-zero. In this situation one is already at maximal twist. However the 
unavoidable 0(a) terms affecting any determination of the critical mass can be 
further tuned in an "optimal way", i.e. in a way such that the residual 0(a 2 ) 
lattice artifacts in physical quantities remain under control as the pion mass is 
decreased. We briefly explain how this can be achieved in practice and to what 
accuracy, following the work of ref. [H]. In the Symanzik expansion of the lattice 
expectation value of a multilocal operator O computed at a bare quark 

mass /J, q there will appear at 0(a 2 ) terms which are proportional to 

eM oc - 2 eM, (4) 

H'q 

where 

^(/i g ) = |(fi|£oddk°(0))|™ nt . (5) 

Here and |7r°(0)) denote the vacuum and the one-pion neutral state at zero 
three-momentum, respectively. With the symbol 

£ 0(W = a£ 5 + a 3 £ 7 + ... (6) 

we indicate the set of operators of odd dimension in the Symanzik local effective 
Lagrangian that describes the maximally twisted lattice theory. From eq. (OD) one 
recognizes that cut-off effects may become large when m 2 gets small. 

The general strategy to avoid these large cut-off effects is to tune £,r(/x 9 ) to 
zero, or at least to reduce it to 0(am 2 AQ CD ) by adjusting the value of K crit . 
One way to realise this is precisely to tune mpcAc — as explained above. In 
particular it is sufficient to impose the vanishing of the PCAC mass at fi q = 
A^.min P3] • An analysis a la Symanzik of the correlator in the numerator of eq. ([3]) 
shows that, if k is such that mpcAc vanishes at a given value of [i q (provided 
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fig < Aq CD ), then ^(/ig) is suppressed in a sufficiently strong way, namely one 
gets (note that £ n has mass dimension three) 

= 0(afi q A 3 QCD ) + 0(afx 2 q A 2 QCD ) + 0(a 3 A 6 Q CD ) , (7) 

and thus (see eq. flU)) 



UN) ^ A N , ^/ N , ^ A ° 2A3 



0(aA QCD ) + 0(a/i 9 ) + 0(aA QCD — ^) • (8) 
AV^-qcd ^q 

In this situation, the ratio ^(/Aj) I /AjAq CD remains small as long as \i q > a 2 Aq CD . 

For each value of \i q in the region cl 2 Aq CD < /i 9 < Aq CD <C a" 1 the value of k 
at which mpcAc vanishes provides a legitimate estimate of K cr it and hence of m cr it. 
Estimates of m crit corresponding to different values of \i q differ by 0{a\x q Aqp^) 
from each other. In particular, working at « C rit (A%min) leads to 0(a 2 ) cutoff effects 
which are at worst of the form a 2 (/i gjmin //ig) 2 and thus perfectly tolerable as long 
as A*g ^ A^min > ^Aqcd- This result can be checked by expanding £ n (iA q ) around 
H>q,min in eq. (jEJ) and using the expression of ^(/x 9 , m in) from eq. (JTj). 

1.3 Numerical precision for tuning to maximal twist 

It remains to be discussed to what numerical precision the condition mpcAc = 
has to be fulfilled. This question is important if one wants to avoid that 
numerical uncertainties jeopardize the tuning procedure. Suppose |ampcAc| = 
ae ^ 0, where ae denotes a small deviation, due to numerical limitations, from 
the condition of vanishing PCAC mass. As a rule of thumb the value of ae can 
be taken as the maximum (in modulus) between the finite statistics central value 
of ampcAc and its (estimated) standard deviation. It then follows by expanding 
£?r = £vr (/A?, e) around e = 



6r(AV e) = £«(p q ) + O (A| CD e) 

« 0(afi q A 3 QCD ) + 0(an 2 q A 2 QCD ) + 0(A 2 QCD e" 



Thus for the relative size of £ n compared to the actual value of the quark mass, 
one gets 

^tf^- = O(oAqcd) + 0(a/i g ) + 0(— ) . (10) 
AV*-qcd A i <? 

A smooth approach to the continuum is, of course, guaranteed when \e/fi q \ is of 
order oAqcd or smaller. In fact, from the form of the dimension five term in the 
Symanzik effective Lagrangian of the twisted mass lattice QCD, it follows that, 
close to maximal twist, aAqcD^ / 'fJ, q \ is the expected order of magnitude of the 
(unwanted) relative 0(a) cutoff effects stemming from violations of the condition 
of vanishing PCAC mass. The requirement |e//i g | < oAqcd thus implies that 
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the relative magnitude of these unwanted cutoff effects satisfies the constraint 
a ^q,CT>\ e l a2 ^-QCD' which is consistent with 0(a) improvement. 

In practice, since we are interested in simulations performed at lattice spacings 
about (or slightly below) O.I fm, where oAqcd ~ 0.1, a value of \e/ ji q \ < 0.1 (and 
ideally decreasing with a) will represent an acceptable level of precision in the 
procedure of determining the critical mass. This condition is not too restrictive 
as in actual simulations it is sufficient that it holds at /ig jiri i n . We also remark 
that in order to check scaling and perform a reliable continuum extrapolation, 
the value of fi q , m i n should be kept roughly fixed in physical units as the lattice 
spacing is decreased. 

Although these theoretical arguments show that we can work in conditions 
such that we are effectively left with only 0(a 2 ) lattice artefacts, numerical com- 
putations are required to check the scaling behaviour and determine the order of 
magnitude of the coefficient multiplying a 2 terms for the observables of interest. 
In this paper, where data at only one value of a are analyzed, we cannot evaluate 
these coefficients. Nevertheless, for the observables we discuss here preliminary 
results from our collaboration presented in ref. [To! [To] indicate that the residual 
cutoff effects are indeed small and consistent with 0(a) improvement. 

2 Computations in the charged meson sector 

In this paper we will be mainly using the twisted quark basis where the fermionic 
action takes the form (CQ). Even though there is no fundamental reason for this 
choice, employing the twisted quark basis makes immediately transparent the way 
several computational methods, which have been invented for, or widely applied 
to, untwisted Wilson fermions, carry over to the case of maximally twisted Wilson 
quarks. Of course, in such an unphysical basis, the two flavour components of 
the fermion field \ = (w, d) T appearing in the action do not coincide with the 
canonical quark fields in the "physical" basis, ip = (%,h ys , d p h ys ) T , rather the 
former are related to the latter by the axial rotation 

x = e'^'X^i^ & « = e-^ /4 « phys , d = e l ^% hys , (11) 

which we write here in the case of maximal twist, uj = ir/2. Since the axial 
transformation above is flavour diagonal, the names of the components (u,d) of 
the twisted basis field x are still appropriate to their flavour content. In spite 
of that, the correct interpretation of gauge invariant composite bare operators in 
the (x, x) basis is obtained only once they are expressed in terms of the physical 
basis bare fields (ip, ip). Examples concerning quark bilinear fields can be found 
in Appendix [A] 

In this context it may be useful to remark that, since parity and isospin are 
no longer exact symmetries (recall however that I3, the third isospin component, 
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is unbroken), a physical basis bare composite operator with given formal par- 
ity and isospin properties can interpolate a hadron with opposite parity and/or 
different isospin. As a consequence in the quantum-mechanical representation 
of the correlators there will be contributions containing matrix elements of a 
physical basis composite operator with given formal parity and isospin between 
the vacuum and a state with opposite parity and/or different isospin, as well 
as between a neutral pion state (which has the same lattice quantum numbers 
as the vacuum) and a state with the same parity and isospin properties as the 
considered operator. Such parity- and/or isospin- violating matrix elements are 
of course of order a. Their occurrence in the quantum-mechanical representation 
of correlators is not in contradiction with the 0(a) improvement of the expec- 
tation values of parity-even, or isospin- invariant, multilocal operators [3]. For 
these specific correlators, indeed, an analysis a la Symanzik shows that each term 
of their quantum-mechanical representation can contain only an even number of 
0(a) factors given by parity- and/or isospin- violating matrix elements 0. 

From the formulae in Appendix |A] it is clear that at maximal twist, uj = n/2, 
the operator dj^u is associated to the ir + meson, in the sense that (dj^uy creates 
the 7r + state from the vacuum. The two-point tt + meson correlator receives 
contributions only from (fermionically) connected diagrams, and after integration 
over fermion fields, it is given by 



where (...) means average over the gauge ensemble, the trace tr[. . . ] is restricted 
to spin and colour indices only, and we denote by G u (0,t) the propagator for a 
M-quark from to t, and correspondingly by Gd the similar propagator for the 
d-quark. Here three-space indices are understood as at this stage we need not 
specify the spatial separation, or equivalently the three-momentum. We can use 
the identity □ Gd(y, z) = ^G u (z, |/) + 75 to relate the connected correlator (|T2l to 
propagators from a common source (at time Xq = 0) through 



Thus only propagators for one flavour at one source point are needed for the 
computation of the charged meson correlator. As we discuss later, it is more 
efficient, however, to evaluate correlation functions from a wider set of sources. 

In the Table below we give the correspondence between bilinear operators 
of the form dTu, where T is an hermitian combination of Dirac 7-matrices, and 
the mesonic state that is associated with each of them (in the limit a — > 0, i.e. 
neglecting 0(a) contamination from states of different parity and isospin). 

lr rhis result essentially follows from the property that, at maximal twist, the order a piece 
of the Symanzik effective Lagrangian, aC^, is odd under parity and the flavour exchange u <-> d. 

2 Here (with a little abuse of notation) by + we mean complex conjugation and transposition 
with respect to spin-colour indices only, while y = (y, yo) and z = (z, zq) are the spacetime 
coordinates. 



C(t) = (tr[G u (0,t)75G d (t,0)75]) 



(12) 



C(t) = (ix[G u (0,t)G u (0,t) + ]). 



(13) 
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Meson 


Operator 


7r± 7T ± , 
4 


dj 5 u, J70M, di^olsu 
di'YfioU, diii^u, d^u 
dijqo^u 
du 



In this Table, Xf labels an isotriplet state with J p = + , for which there is no 
experimental candidate. We note that the associated operator is in the continuum 
a component of a conserved current in the theory with two mass degenerate 
quarks. 

We evaluate the two-point (connected) correlators for all the pairs of operators 
in the same line of the Table above. In view of the symmetries of the lattice theory 
at maximal twist [21 E], such correlators are in general non zero: e.g. the correlator 
obtained from the insertion of the first (or second) operator in the second line of 
the Table with the third operator in the same line is an 0(a) quantity (in fact 
p ± and af carry different continuum quantum numbers). Since we also use a 
local and extended (fuzzed) source and sink in all cases we consider, we will have 
either 6 x 6 or 2 x 2 matrices of correlators available. 

Therefore, we measure in general correlation functions of several different 
pairs of operators ((0 Q 0/?), with a, ft = 1, . . . N) at source and sink. We then 
use a factorizing fit expression where % = 1, . . . M states (with energy denoted by 
Ei) are included 

M 

C aP (t) = <44( e "** ± e- E ^ T - t] ) • (14) 
i=i 

Here T is the lattice temporal extent and the ± sign is determined by the prop- 
erties of the chosen operators under time-reflection. By simultaneously fitting 
N x N correlators with M states, we can optimally determine energies and cou- 
plings. From them we evaluate other quantities of interest, such as af n and 
ompcAc- We use conventional methods to determine the optimal t range, N- and 
M-values to be employed in the fits. We take into account statistical correlations 
among observables [T7j through correlated fits to establish that the x 2 value is 
acceptable. Our final fitted values are obtained from uncorrelated fits, since that 
introduces less bias p2], although the x 2 values are smaller than those obtained 
including correlations. We also checked that the fits are stable when taking into 
account correlations. For pseudoscalar mesons we use mainly M = 1 as well as 
N = 4 or 6, and select the minimum value of t such that the effective masses (or 
energies) from different matrix elements agree. 

We conclude by recalling that, owing to reflection invariance of the lattice 
action (the Euclidean analog of the Minkowski complex-conjugation, a symmetry 
that is preserved by Wilson fermions, either chirally twisted or not, see ref. [3]), 
all the correlators that are expectation values of fields with definite reflection 
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properties are either real or purely imaginary, depending on whether the whole 
field product has even or odd reflection-parity. In particular, the expectation 
values of multilocal fields with negative spatial parity, which are 0(a) quantities, 
come out to be purely imaginary if one does not take care of inserting the i-factors 
that are needed to render the multilocal field even (rather than odd) under the 
reflection. 

2.1 Quark propagators from stochastic sources 

Although it is feasible to use w-quark propagators from 12 colour-spin sources 
(with each source being non-zero only for one colour-spin combination) at one 
space-time point to evaluate mesonic correlators, it is preferable to use the in- 
formation contained in the gauge configurations more fully, especially in the case 
of such CPU-time expensive simulations as dynamical quark simulations. One 
efficient way to achieve this goal is to use stochastic sources. To keep the noise-to- 
signal ratio reasonable, it is mandatory to use time-slice sources rather than full 
volume sources. A great reduction of the noise-to-signal ratio over conventional 
stochastic methods (see ref. [TS] for a review) can be obtained [T9"} |2"U] by using 
the "one-end-trick" which is described below. A similar method, called random 
wall, was used by MILC pi]. 

The starting point of all stochastic methods for computing quark propagators 
is the introduction of random sources, £[, where i = 1, . . . V s spans the set of the 
source degrees of freedom (colour, spin, space, time) and r = 1, . . . R labels the 
noise samples generated for each gauge configuration. The corresponding average 
satisfies 

lim [£?£,■]* = Ay, lim = , (15) 

K— >oo H — >oo 

which can be achieved by various different noise choices, such as = (±l±i) / \[2 
or gaussian (complex) noise. 

As a next step, we invert the lattice Dirac matrix M (for one given quark 
flavour) on each sample of this source, 

% = MTk?k, ( 16 ) 
so that averaging over r (R samples) gives 

= I^M^Qr = MJ, 1 + noise , (17) 

where j can be arbitrary and i belongs to the set of indices for which the source 
is no n- vanishing, which we assume to be of size V s . The quantity (JTTJ) is an 
unbiased estimator of the quark propagator from i to j. Unfortunately, here the 
noise is expected to be as ~ \fV s j VR whereas the signal is ~ 1 at best. Variance 
reduction is thus very necessary. Furthermore for a meson correlator, the signal 
behaves as exp(— m mcson t) which decreases rapidly with increasing t. 
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The 'one-end-trick' allows P3JJ [20] a more favourable signal-to-noise ratio. 
Consider the product where the stochastic source is now non-zero for all 

colour-spin indices and all space points at only one time, denoted by to (time-slice 
source). Then upon averaging over r one has 

[WjU = [(M^aTM^tlU = M„ U M]? + noise, (18) 

where the sum over k includes all source components. This quantity is an unbiased 
estimator for the product of the quark propagators Mj~ k M ki + from the source to 
sites % and j on each gauge configuration. Then contracting with 8ij and summing 
over space at fixed time-slice t yields the full zero three-momentum (^-channel) 
correlator from to to t. The noise counting is now more favourable. There are 
noise terms, which yield a standard deviation of order V s , but the signal itself is 
of order V s . This is such big an advantage that it is sufficient to employ just one 
sample of noise per gauge configuration (R = 1). As we discuss below, the optimal 
way to choose the time-slice (to) at which the stochastic source is located, is to 
change it randomly as the gauge configuration is changed. It should be remarked 
that the 'one-end-trick', as formulated above, only works for the case of a zero 
three-momentum interpolating field of the form d'-f^u at the source time (to). 

A convenient extension of the 'one-end-trick', that allows meson-to-meson 
correlators with any Dirac structure at the source to be evaluated, requires con- 
sideration of four {(3 = 1, 2, 3, 4) "linked" sources of the form 

t(/3;to) x x n 

where a and c are Dirac and colour indices respectively, while 77 is a non-vanishing 
noise field. Such sources, which are non-zero only on a given time-slice (to) and 
when the Dirac index value equals /3, are called "linked" because they involve a 
common noise field rj. [f| One can check that by replacing £ and £* in the l.h.s. 
of eq. (|18|) by two of these linked sources, say £^ ; '°) and £vy>'*°)* : and choosing 
appropriately (3 and 7, it is possible to evaluate the two-point correlators with a 
field of the form dTu at the source (xq = to) with any Dirac matrix T. This very 
useful extension, which we have thoroughly exploited in the present paper, comes 
at a moderate price. One must in fact only perform four separate inversions (per 
gauge configuration and per noise sample), one for each of the four linked sources 

£(/3;to) ; /3=1 ) ...4. 

To further extend the one-end trick with linked sources to non-zero three- 
momentum or to spatially non-local mesonic operators is completely straight- 
forward, at the cost of more inversions. One creates further linked sources Ft; 
(where F denotes a product of links) with the desired spatial properties, and 
computes the quark propagators originating from them, <fip = M _1 F£. Combin- 
ing the latter with the quark propagator stemming from £, i.e. <fi = M -1 £, yields 

3 Note that "linked" sources are different than "spin-diluted" sources [HI [22] since these 
require different random numbers for each spin. 
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the product (jf(j)Fi from which, upon averaging over the noise, one can evaluate 
a set of correlators with the meson field d'j^Fu inserted at the source (and all 
possible spatial structures at the sink). Employing linked sources, as explained 
above, one can finally evaluate correlators with the meson field dTFu inserted at 
the source with any spatial structure F and Dirac matrix T, while retaining all 
advantages of the one-end trick. 

In this work we use fuzzing, see Appendix ID . 2 1 and ref. [23], to create spatially 
non-local meson operators, since this procedure is computationally fast also at 
the sink. The fuzzed meson source is constructed from a sum of straight paths 
of length 6a, in the six spatial directions, between quark and anti-quark. These 
straight paths are products of fuzzed gauge links. Here for the fuzzed links we 
use the iterative procedure defined in Appendix ID. 21 with X s = 0.25 and n = 5. 

In principle one could hope to extend the approach described above to bary- 
onic correlators (choosing £ as a cubic root of 1) but the signal to noise ratio will 
be less favourable (noise induced standard deviation will be ~ V s versus signal 
~ V s ). Unfortunately one finds that this extension of the stochastic method to 
baryons is not any improvement over using point-like sources. In general, the 
choice of the optimal stochastic methods needs to be investigated on a case by 
case basis. 

2.2 On the way of choosing the source time-slice 

As discussed above, we invert on spatial-volume stochastic sources located at 
time to, where < to < T can be chosen differently for each gauge configuration. 
We have explored two ways of changing the source time-slice to- One consists in 
moving to cyclically through the lattice. This means that we choose n equally 
spaced values for the source time locations, t , < i < n. Then we invert on 
the j'-th gauge configuration using sources that are non-vanishing only at the 
time-slices to — t$ mod . Hence, we invert from the same time-slice only every n 
configurations, i.e. after one cycle. Even though this method should decorrelate 
the measurement on two consecutive gauge configurations better than when the 
time-slices are kept fixed, it has the drawback that after a relatively short number 
of configurations the same time-slice is used again. Actually, at least for the 
mesonic correlators studied in this paper, it turns out that two measurements 
from the same time-slice, but 8 trajectories apart, are much more correlated than 
two measurements from different time-slices, but only two trajectories apart. 
Furthermore the analysis with the V method of ref. [21] described in sect. 14. II and 
Appendix O becomes ill-defined, because translational invariance is broken. This 
invariance can be recovered, however, by averaging over cycles and using the V 
method on the cycle-averaged ensemble. 

The second way of moving the time-slice we explored was to choose the value 
of to randomly for every gauge configuration we inverted on. This method also 
maintains translational invariance properly for a large enough configuration en- 
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semble. It is therefore expected to work better than the aforementioned cyclical 
way. This will indeed turn out to be the case, as we shall see below, where we 
discuss in more detail the effects of these two ways of generating source time- 
slices. 

3 Computations in the neutral meson sector 

Lattice QCD with maximally twisted Wilson fermions enjoys the remarkable 
property that, even if the action is not 0(a) improved, all the physically rel- 
evant observables are affected by cutoff effects only at order a 2 (and higher). 
Among these 0(a 2 ) cutoff effects will be a violation of parity and (in part) 
isospin. Isospin and parity violations have several consequences for meson spec- 
troscopy. For instance 1) neutral and charged mesons can have different masses, 

2) quark-disconnected contributions are needed for neutral isovector mesons and 

3) correlators receive contributions from states that in the continuum limit carry 
different parity and isospin quantum numbers [3J. Here we discuss how we com- 
pute the correlators for neutral mesons and, in particular, the quark-disconnected 
(for brevity called simply "disconnected" below) contributions. We illustrate our 
approach in the relevant case of the neutral pseudo-scalar meson. 

The neutral pion can be created by the operator y^ip'j^Ts^p which, at maximal 
twist, in the twisted quark basis reads (i/VtyxX = {i/V2){uu + dd). When this 
operator is inserted at source and sink, we will have to consider the correlators 

G tot (t) = ((uu + dd)(t)(uu + Jd)(0))/2 , (19) 

where again three-space indices are understood. The latter can be rewritten in 
the form 

C tot (t) = C(t) + D(t) , (20) 
C(t) = -(tr[G u (0,t)G tt (t,0)] + tr[G d (0,t)G d (t,0)]>/2, (21) 
D(t) = (tr[G u (0, 0) + G d (0, 0)]tr[(G u (t, t) + G d (t, t)])/2 , (22) 

with the trace tr[. . . ] running only over spin and colour indices. As usual, we can 
relate the connected contribution (C) to propagators from a common source (at 
time Xq = 0) through 

C(t) = -(tr[ 75 G M (0,t)7 5 G d (0,t) + ] +tr[ 75 G d (0,t)7 5 G u (0,t) + ])/2. (23) 

The disconnected contribution can be expressed as 

D(t) = (tr[G u (0, 0) + G u (0, 0) + ]tr[G u (t, t) + G u (t, t) + }}/2 . (24) 

Thus we see that to evaluate the correlation ( fl9j) we need both u and d-quark 
sources for the connected contribution as well as an evaluation of the disconnected 
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contribution for w-quarks at both initial and final t-value. This is at variance 
with the 7r + correlator which can be evaluated from a w-quark source alone and 
which has no disconnected contribution. The evaluation of the disconnected loops 
is detailed in Appendix [Bj including discussion of both the hopping-parameter 
method for the reduction of the stochastic noise [25] and a new powerful method 
of variance reduction applicable in many cases. 

In the table below we give the correspondence between bilinear operators of 
the form uTu ± dTd, where T is an hermitian combination of 7-matrices, and the 
neutral mesonic state that is associated with each of them in the limit a — > (i.e. 
ignoring 0(a) contaminations from states of different parity and isospin). 



Meson 


Operator 


A A /o 
77, 77, a° 

uj, uj, b\ 
a? 
fi 

xl 

X° 


Xilol5T 3 x, XX, Xl^zX 
XH0I5X, XT3X, X15X 
XliT3X, XiHlol5X, XhiloT3X 

xiiX, xhao^nx, xhaox 

Xililhnx 
XililbX 

X10T3X 

XloX 



Here X® (Xq) labels an isotriplet (isosinglet) state with J PC = + ~, for 
which no experimental candidate is known. We remark that these operators are 
conserved isotriplet (isosinglet) currents in the continuum theory with two mass 
degenerate quarks. 

As in the charged channel, we evaluate the two-point correlators where only 
pairs of meson operators appearing in the same line of the Table above are in- 
serted. Since we use a local and extended (fuzzed) source and sink in each case, 
we have either 6 x 6 or 2 x 2 matrices of correlators available. The connected 
correlators are actually the same for certain states of different isospin (e.g. rj or 
7r). The same 'one-end-trick' discussed above, based on the use of stochastic 
time-slice sources with random choice of its position on each gauge configuration, 
can be used for the connected neutral correlator. 

In more detail, we use four "linked" sources see sect. 12. ip and further 

four fuzzed sources based on the same noise field (F£^) to compute ordinary 
and fuzzed w-quark propagators from one time-slice to all points. This set of 
eight sources is just the same we used to evaluate correlators of charged mesons. 
For neutral mesons, we inverted the lattice Dirac matrix of the d-quark on the 
same (non-fuzzed) four stochastic sources (£^) as above and on the correspond- 
ing four stochastic sources with the lowest possible three-momentum (2tt/L, for 
simplicity taken always in the x-direction) . In principle mesonic operators with 
non-zero anisotropic three-momentum have less symmetry then their counter- 
parts with vanishing three-momentum, implying that more correlators may take 
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Run 


L 3 xT 


a/lq 


iVtraj 


iVcfg 


Bla 


24 3 x 48 


0.0040 


5000 


2500 


B±b 




0.0040 


1341 


670 


Bic 




0.0040 


3380 


1690 


B 2 




0.0064 


5192 


2500 


B^a 




0.0085 


3753 


1876 


Bm 




0.0085 


940 


470 


Bi 




0.0100 


5000 


2500 


B5a,b 




0.0150 


2500 


1250 



Table 1: Summary of all simulation points. We give the lattice size L 3 x T and the value of 
the twisted mass a[i q . In the last two columns we quote the number of equilibrated trajectories 
-^traj produced and the number of configurations A^fg saved to disk and finally stored within 
ILDG, see the review [27] for further links and references. All runs listed in the Table have 
been performed at (3 = 3.9 and k ~ 0.160856. 

non-zero values. Here we do not evaluate these additional correlators. We do take 
care, however, to distinguish between the vector meson correlators with three- 
momentum parallel to spin and those with three-momentum perpendicular to it. 
As shown elsewhere [2E], a study of the difference between these correlators can 
shed some light on the mixing of p mesons with their decay products {itix). 

4 Simulation algorithm and error analysis 

In this section we provide details on the algorithms we used to generate the 
gauge configurations and information on the methods employed for the estimate 
of statistical errors. 

In Tabled] we give the list of the key parameters characterising the simulations 
we are going to use in this paper. All simulations Bi — B 5 have been performed 
at a fixed value of the gauge coupling (3 = 3.9 and a fixed value of the hopping 
parameter k = (8 + 2amo)~ l = 0.160856 on 24 3 x 48 lattices. In addition to the 
values of ap q we provide in Table [1] the number of trajectories, A^raj, produced 
after allowing for 1500 equilibration trajectories, and the number of gauge config- 
urations, iV c f g , that were saved on disk (one every second trajectory). For every 
value of a/i g we have reached ~ 5000 equilibrated trajectories. 

In case we have several ensembles (as for instance for B{) or several replicas (as 
for instance for B§) for the same lattice parameter set, we denote this by adding 
an extra subscript, a, b, .... For our smallest value a\i q = 0.004 we extended our 
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statistic from about 5000 trajectories (ensemble B la ) to ~ 10000 trajectories (if 
also trajectories from ensembles Bi b and Bi c are counted). 

The algorithm we used is a HMC algorithm [28] with mass preconditioning [29l 
[30] and multiple time scale integration, as described in detail in refs. [3U [32] • The 
algorithm parameters we employed for the various runs can be found in Table [2} 
where we mostly follow the notation of ref. [21] • The integration schemes we used 
are the Sexton- Weingarten (SW) scheme [33], the second order minimal norm 
scheme (2MN) [M] and its position version (2MNp). We also list the number of 
integration steps Ni for time-scale i (for details see ref. [H]). We recall that N 2 
represents the number of integration steps of the outermost (largest) time-scale. 
Thus the number of integration steps of the smallest (i.e. innermost) time-scale 
(the one referring to the gauge field integration) is given by N 2 ■ N\ ■ N . The 
preconditioning mass is given by jli = 2k/j,i, with fix typically larger than \x q by 
a factor O(10). 

The second order minimal norm integration scheme on time-scale i is 
parametrised by one real number, \. We also give the number, iV J csg , of so- 
lutions of the Dirac equation we save for the chronological solver guess [33] with 
the purpose of evaluating the two force terms (i = 1,2) associated to pseudo- 
fermion integration (i = refers to the pure gauge force). The notation iVff = 
means that no chronological solver guess was used there. Finally, we quote the 
acceptance rate P acc observed in the simulation and the integrated autocorrela- 
tion time r int (P) of the plaquette expectation value. The trajectory length was 
set to r = 1/2 in all our runs and we always used Npp = 2 pseudo-fermion fields. 
For details on the linear solvers we employed to invert the Dirac matrix we refer 
to ref. [36] . 

To give guidance on the computational cost of such simulations, we specify the 
resource used at our lightest /i^-value where the CG iterations for one trajectory 
cost about 115 Tflop. The production of 5000 trajectories amounted to about 17 
rack days on the BlueGene/L installation in Julich, with our code running with 
an efficiency of about 18% for the B\ parameter set. 

4.1 Statistical error analysis 

A reliable estimate of the statistical errors on the measured quantities is extremely 
important for many reasons. We discuss here only the points which are of special 
relevance in our analysis. If the basic systematic effects in the lattice simulation, 
originating from the lattice discretization, the finite volume and the mass of the 
dynamical quarks are to be addressed, the statistical accuracy on all the relevant 
quantities has to be understood very well. In fact, on the one hand, relevant 
but tiny systematic effects can only be detected with high statistical accuracy, 
on the other hand underestimated statistical errors can artificially increase the 
significance of systematic effects. The PCAC quark mass, though not a physical 
quantity, plays here a special role, since the precision by which it is set to zero 
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Run 


Int. 


#0,1,2 


fa 


Ao,l,2 


/V csg 

iv l,2 


p 

1 acc 


7int(P) 


B\a,b 


2MNp 


2,3,6 


0.018 


0.19,0.20,0.21 


0,0 


0.85 


47(15) 




SW 


2,3,6 


0.018 


— 


0,0 


0.90 


43(15) 


B 2 


2MNp 


2,3,6 


0.025 


0.19,0.20,0.21 


0,0 


0.90 


23(7) 


Bza,b 


2MN 


2,3,5 


0.020 


0.19,0.20,0.21 


7,1 


0.90 


13(3) 


B 4 


2MNp 


2,3,6 


0.035 


0.19,0.20,0.21 


0,0 


0.90 


15(4) 


Bz, a ,b 


SW 


2,2,6 


0.050 




0,0 


0.90 


30(8) 



Table 2: HMC algorithm parameters. For all ensembles we specify the integration scheme, the 
number of time steps on each time scale No,i,2, the precondition mass fix = 2k/j,i, the A- values 
for the 2MN integration scheme, the number of saved solutions N csg for the chronological solver 
guess, the acceptance rate P acc observed in the run and the integrated autocorrelation time of 
the plaquette r int (P). The trajectory length was set to t = 1/2 for all runs and we used always 
-ZVpp = 2 pseudo-fermion fields. 

is related to the accuracy (see sects. ll.2Hl.3ft by which we can expect to be at 
maximal twist. Secondly, small statistical errors in low-energy hadronic quantities 
is an expected virtue of the twisted mass formulation, where a sharp infra-red 
cut-off ensures a stable MC evolution of the lattice system. Of course we have 
to make sure that an apparently small statistical error does not come as a result 
of large unnoticed autocorrelations in the MC history. So autocorrelations in the 
measured quantities must be accurately analysed. Finally, a detailed analysis of 
the statistical errors delivers as a by-product the integrated autocorrelation time 
ri n t of the studied observable, from which the efficiency of the employed algorithm 
as a function of the simulation parameters can be quantified (see sect. 14. 2D . 

Given the importance of getting a reliable estimate of statistical errors, results 
have been cross-checked using different approaches. As for the estimate of auto- 
correlation times two different kinds of analyses have been performed: one based 
on a standard data-blocking (or binning) procedure and another one relying on 
the so-called T-method [371 El]. m order to keep self-contained this paper we dis- 
cuss these methods in some detail in Appendix O Since there are arguments [21] 
supporting the superiority of the T-method over data-blocking, the former will 
be our method of choice in the evaluation of r int and the error on it, for several 
observables. In particular the T-method has been used to estimate the statistical 
error on the plaquette and mpcAc, which turn out to have large autocorrelation 
times. For all the other observables having significantly smaller autocorrelation 
times data-blocking and T-method typically give quite similar error estimates. 

Cross-correlations among different observables are properly taken into account 
in our error analysis by using standard jackknife or bootstrap [38] or performing 
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fits based on a definition of \ 2 that involves the inverse covariance matrix (see 
eq. (H2l) and the discussion of Method A in sect. 17.11) . 

4.2 Autocorrelation times 

For a primary observable O, i.e. one that can be viewed as a linear combination of 
expectation values of multi-local operators, the integrated autocorrelation time 
is in principle given by 



where T (n) is the autocorrelation function of the observable O (see eq. (EH]))- 
The autocorrelation times for the plaquette and fermionic quantities, like amps, 
a/ps and ampcAc, were determined using the T-method as described in ref. [21] 
(see also Appendix [C]) . This method allows the determination of T- mt also for 
non-primary quantities, as the aforementioned fermionic observables. The values 
for the plaquette integrated autocorrelation time are collected in Table [21 those 
for amps, a/ps and ampcAc i n Table [31 All quoted values are given in units of 
trajectories of length 1/2. 

In the case of the ensemble B\ a we employed the two ways of moving the 
stochastic source through the lattice described in sect. 12.21 As can be seen in 
Table El indeed the random way performs better. This is especially significant 
for ampcACi for which we observe the longest autocorrelation time among the 
fermionic quantities. For amps and a/ps the difference between the two methods 
is not significant. The somewhat larger autocorrelation time of B^ c in particular 
for mps and fp$, stems presumably from the fact that the time slice sources were 
chosen closer to each other than at the other ensembles. 

Table [H gives details on the computational methods employed to extract the 
various fermionic quantities. In the case where the random way of moving the 
source is used, the value of t p reported there represents the number of trajectories 
(of length 1/2) between two consecutively measured gauge configurations. We 
also give in this Table the value of x 2 /d.o.f. obtained when the charged pion 
correlators (for Euclidean time separations in the range 10 < t/a < 23) are fitted 
using the Ansatz (TH|) . 

Looking at the three fermionic observables reported in Table El we observe 
that the integrated autocorrelation times of ampcAc are significantly larger than 
those of amps and a/ps- We attribute the large value of the autocorrelation 
time of ampcAc to the peculiar phase structure of twisted mass lattice QCD with 
Wilson type quarks as discussed in ref. [7J . The simulated values of \i q are not in 
a region where the phase transition occurs. However, the system may still feel the 
presence of this phase transition. The situation is similar for the plaquette value, 
as also discussed in ref. [7], and indeed the integrated autocorrelation times for 
the plaquette and ampcAc are rather similar. 




(25) 
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Figure 1: The Monte Carlo history of the ratio of correlators defining the PCAC quark mass 
estimator described in the text on the configurations of the ensemble B±. The configuration 
number corresponds to the number of trajectories divided by two. 

We show the Monte Carlo (MC) history of an estimator of ompcAc, see 
eq. (1321) . for our lightest quark mass (a/i q = 0.0040) in fig. [TJ More pre- 
cisely, we plot for each gauge configuration the axial-pseudoscalar correlator at 
t/a = 10 (where the pion ground state is dominant) multiplied by the factor 
0.5amps/Cpp(10), where the average over all gauge configurations is used in 
Cpp(10) 0. The plot in fig. [T] shows long-ranged fluctuations in MC time. The 
autocorrelation function in eq. ( 1691) from this data set is reported in fig. [2j 

From it an integrated autocorrelation equal to 32(9) trajectories is obtained, 
in agreement with the result quoted in Table. [H] (third line). 

Concerning the neutral pseudoscalar meson, a study of the corresponding cor- 
relators indicates that the autocorrelations are definitely shorter than 100 (length 
1/2) trajectories. Thus our error estimates, coming from a bootstrap analysis on 
blocked data with blocks made of measurements taken from 80 trajectories, are 
expected to be reliable. 

4 Note that the average of this quantity over all gauges is not our best estimator of ompcAC, 
since it does not exploit the possibility of averaging over t (the Euclidean time separation) . 
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Figure 2: The autocorrelation function (see eq. from the data presented in fig. [TJ The 

vertical line shows the window W from eq. (I70|) used to evaluate the integrated autocorrelation 
function. 



Within the relatively large errors of our estimates of the autocorrelation times, 
it is actually not possible to find a significant dependence on the value of the 
twisted mass a\x q for any of the fermionic quantities discussed here. 

5 The scale from the static potential 

A convenient way to set the scale in lattice simulations is through measurements 
of the static potential and the associated hadronic scale r$ [39] . Although we will 
finally not use ro to set the scale in our dynamical simulations, we will use it for 
a scaling analysis towards the continuum limit and its reliable determination is 
therefore important for us. The scale r is defined via the force between static 
quarks at intermediate distance 

r 2 F(r ) = 1.65, (26) 

where numerical calculations are most reliable and hence are expected to lead 
to very accurate results. We measure the static quark-antiquark potential by 
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Run 



Tmt(cOTlps) 



Tint(a/ps) 



Tint (ampcAc) 



r>cyc 
D \a 


7(1) 


13(4) 


60(24) 


Rrnd 


6.6(1.1) 


8(1) 


20(5) 


prnd 
D \a,b,c 


5.9(7) 


7(1) 


23(5) 


E cyc 


17(4) 


33(8) 


43(14) 


ocyc 
n 3 


10(2) 


11(2) 


66(27) 


R cyc 


7(2) 


14(4) 


54(23) 


-°5a,6 


20(6) 


14(3) 


105(51) 



Table 3: Estimated integrated autocorrelation times for amps, a/ps & n d awipcAC- The labels 
eye and rnd refer to the cyclic and random choice of the source, see text. All integrated 
autocorrelation times are given in units of trajectories of length 1/2. The fact that for the 
ensemble B§ yc we find a rather large autocorrelation time with, however, a large error we 
attribute to the usage of 2 replica in the analysis. 



determining expectation values of Wilson loops of size r x t on our ensembles of 
configurations. Unfortunately, the relative errors of the Wilson loop expectation 
values increase exponentially with the temporal extension t. To reduce these 
statistical fluctuations one can employ improved static actions amounting to use 
modified temporal links for building Wilson loops[E However, it is also important 
to enhance the overlap with the physical ground state of the static system and 
this can be achieved by invoking iterative spatial smearing techniques together 
with a variational method to extract the ground state. The computational details 
for calculating the static potential are given in Appendix [D] while in the following 
we want to concentrate on analysis details and physical results. 



5.1 Analysis details and results 

In order to extract the physical scale through eq. ( 1261) we need an interpolation of 
the potential and correspondingly of the force between the quarks for arbitrary 
distances r. This interpolation is achieved by fitting the form of V(r) with the 
ansatz M 

a 

V(r) = V + - + ar. (27) 
r 

We employ a two step procedure to perform the interpolation. First we extract 
the values of the potential V(r) for each r separately using standard variational 
techniques. In a second step we fit directly the potential ansatz in eq. (1271) to 

5 See ref. [ID] for a first use of this idea. 

6 Note that we do not use tree level improved distances. 
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Run 


method 


t P 


timeslices 


X 2 /d.o.f. 




cyclic 


8 


0, 12,24,36 


0.12/39 


B\ a ,b,c 


random 


10 


- 


2.50/39 


B 2 


cyclic 


16 


0,6,12,18,24,30,36,42 


1.15/39 


B 3 


cyclic 


8 


0,12,24,36 


2.38/39 


B 4 


cyclic 


8 


0,12,24,36 


1.85/39 


B<j a ,b 


cyclic 


8 


0,12,24,36 


0.99/39 



Table 4: Measurement methods for fermionic quantities. The source timeslices are either 
chosen in a cyclic way or randomly. For the cyclic way, t p denotes the number of trajectories 
between two configurations for which the same time-slice was used. For the random way, it 
specifies the number of trajectories between two measured configurations. In the cyclic case 
we also specify the time-slices where the source was located. Finally, the fit range we chose 
to determine a/ps and amps was always 10-23 and the corresponding values of x 2 /d.o.f. for a 
2x2 factorising fit (sec cq. (fT4")0 are quoted in the last column. 

the Wilson loop correlators taking into account all spatial and temporal cross- 
correlations in the data. These two steps are now described in more detail (see 
also [Hj). 

We use five spatial smearing levels S n U, n = 8, 16, 24, 32, 40, and hence we 
end up measuring a 5 x 5 correlation matrix C (r, t) of spatially smeared and 
temporally improved Wilson loops (see Appendix 0) . The variational method 
results in a linear combination of the string operators, which projects sufficiently 
well to the ground state of the string, i.e. has the effect of eliminating the closest 
excited string states. This is done by solving the generalised eigenvalue problem 

C(r,t 1 )v i = \ i (r-,t ,t 1 )C(r,t )vi, A, > . . . > A 5 , (28) 

with to = 3a and t\ = 4a and projecting the correlation matrix to the eigenspace 
corresponding to the largest eigenvalue, i.e. the ground state, 

C(r,t) = (v 1 ,C(r,t)v 1 ). (29) 

Based on effective masses and on a x 2_ t es t which takes the temporal cross- 
correlations between C(r,t) and C(r,t') into account, we choose a plateau region 
from £ min to t m ax- Too small t values distort the results due to contamination of 
excited states, while too large values introduce noise. Examples of effective mass 
plateaus and the chosen fit ranges are provided in figure [3] for the Wilson loop 
correlators of the five ensembles -E>i,...,5 at quark-antiquark separation r/a = 4. 

The results of our fits are collected in Tables [5H9] where we list the plateau 
regions (fit range), the values of the extracted potential, V(r), and the x 2 P er 
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Figure 3: Effective masses for the ground state energy of the Wilson loop at quark- antiquark 
separation r/a — 4, for ensembles B\ to B$. 



degree of freedom. The uncertainties in the extracted values of V(r) are calculated 
using a non-parametric bootstrap method [38]. 

This first step allows to determine for each r the value of the potential V(r). 
A straightforward strategy to evaluate the scale r^/a is then to fit the numbers 



22 



r/a 


^min 


^max 


aV{r) 


X 2 /d.o.f. 


4 


9 


11 


0.3308(10) 


1.09 


5 


8 


10 


0.3974(15) 


0.09 


6 


7 


9 


4577(26) 


2.48 


7 


7 


9 


0.5126(31) 


1.60 




Table 5: Fit parameters for the Wilson 


loop correlators for run B\. 




r/a 


^min 


^max 


aV{r) 


xVd.o.f. 


4 


9 


11 


0.3319(07) 


1.56 


5 


8 


10 


0.3994(09) 


3.27 


6 


7 


9 


0.4589(12) 


0.16 


7 


7 


9 


0.5129(21) 


0.13 



Table 6: Fit parameters for the Wilson loop correlators for run B^. 



so obtained with the ansatz ( ]27l) . However, one can diminish the errors on the fit 
parameters by exploiting the fact that data at different values of r are correlated. 
Therefore, in a second step we use the ground state projected correlator C(r,t) 
to estimate the covariance matrix 

Cov(r,t;r',0 = (C(r, t)C(r>, f)> - (C(r, t))(C(r', t')) 

from the bootstrap samples of C(r,t) and use Cov(r, t; r', t') to construct the \ 2 
function (see the discussion in sect. 17. ip . The r and t dependence of C(r,t) is 
fitted with the formula (see eq. 1 )2*71) ) 

C(r, t) ~ Z(r) exp [-tV(r)\ = Z(r) exp [-t (V + a/r + or)] . (30) 

For the temporal fit interval we use the fit ranges t m i n (r) to t max (r) determined 
in the first step. The fit range in r is chosen so as to include only a few values of 
r closest to ro in order to minimise both the statistical error and the systematic 
error coming from the choice of the interpolation formula. Once the best fit 
parameters (V , a and a) in eq. ff3H|) are found, the value of r /a is obtained 
straightforwardly by computing the static force from the derivative w.r.t. r of 
eq. ( 1271) and imposing the condition f)2~6"]) that defines r . 

A compilation of the results of our fits is provided in Table [1(3] where we give 
the number of measurements N meas and the x 2 P er degree of freedom, in addition 
to the results for r /a. The final error on r Q /a is estimated through jackknife 
and bootstrap procedures using binning to take residual autocorrelations into 
account. In Table [10] we give the errors from the jackknife procedure using a 
binning factor equal to 4. 
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r/a 


^min 


^max 


aV{r) 


X 2 /d.o.f. 


4 


9 


11 


0.3321(08) 


0.75 


5 


8 


10 


0.4002(11) 


4.39 


6 


7 


9 


4617(17) 


0.40 


7 


7 


9 


0.5177(22) 


0.02 




Table 7: Fit parameters for the Wilson 


loop correlators for run B3. 




r/a 


^min 


^max 


aV(r) 


xVd.o.f. 


4 


9 


11 


0.3335(06) 


1.40 


5 


8 


10 


0.4013(12) 


0.41 


6 


7 


9 


0.4616(15) 


0.37 


7 


7 


9 


0.5177(25) 


0.98 



Table 8: Fit parameters for the Wilson loop correlators for run B4. 



5.2 Discussion 

There are several sources of systematic effects which can distort a precise and 
accurate determination of the scale r /a. Here we would like to discuss a few 
checks that we have performed in order to asses these systematic effects and 
some procedures to minimise their influence. 

Excited states 

First of all there are contaminations of the ground state energy of the Wilson 
loops from excited states. We expect that these should be eliminated by our 
variational calculation of the ground state and our choice of the fit range in t, 
and we have carefully checked the stability of the results under variation of the 
fit parameters (see figure H]). In particular we have checked that we can resolve 
the first excited state and that the ground state energy remains stable under this 
procedure. Moreover we have also checked the stability of the ground state under 
a truncation of the variational operator basis. We would also like to point out 
that the fit ranges in t were not chosen independently for each value of \i q and r, 
rather we chose them after taking a global view of the effective mass data for all 
values of fi q at given fixed r (see fig. [3] for the case r = 4a). This procedure makes 
sense since the /^-dependence of the Wilson loop correlators is expected to be 
rather weak (see below) and is particularly useful in cases where the choice of 
the fit range for the effective masses cannot be determined unambiguously given 
the available statistics. Finally we note that contaminations from excited states 
tend to increase the potential energies and the effect will be more pronounced for 
the larger Wilson loops. As a consequence, residual contributions from excited 
states will tend to decrease the value of r /a. 
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r/a 


train 


^max 


aV{r) 


xVd.o.f. 


4 


9 


11 


0.3373(07) 


1.01 


5 


8 


10 


0.4062(08) 


3.21 


6 


7 


9 


4692(14) 


1.24 


7 


7 


9 


0.5272(22) 


2.28 




Table 9: Fit parameters for the Wilson 


loop correlators for run Bg, 




Run 


-^meas 


d.o.f. 


X 2 /d.o.f. 


r /a 


B 1 


625 


5 


1.44 


5.196(28) 


B 2 


695 


5 


1.80 


5.216(27) 


B 3 


598 


5 


2.58 


5.130(28) 


B 4 


602 


5 


0.57 


5.143(25) 


B 5 


645 


5 


1.92 


5.038(24) 



Table 10: Results of the fits for the scale ro/a from the static potential. The fit range was 
always r/a = 4 — 7. The number of measurements, -ZV moas , and % 2 /d.o.f. are also reported. 



Interpolation error 

The interpolation of the potential (or the force) as a function of r is not unique. 
Here we would like to emphasise that we use eq. (12?]) only locally as a simple 
interpolation ansatz without attaching to it any special physical meaning. As 
a check of this interpolation ansatz, one can use separately the matrices of 
correlators computed for r/a = 4 — 6 and for 5 — 7 to obtain two different 
determinations of ro/a. Their difference then provides an estimate of the error 
coming from the interpolation procedure. It turns out that our choice of the fit 
range r/a = 4 — 7 covers this spread typically within 1-2 standard deviations of 
our final result (see figure HJ) . 

Correlations 

We have already pointed out that it is important to take both the spatial 
and temporal cross- correlations of the Wilson loop operators into account 
when fitting them to the ansatz (1301) . Our finite statistics limits ourselves 
to short fit ranges in order to obtain a stable covariance matrix, and this is 
one of the motivations for the rather narrow fit ranges in t in Tables EH9] 
In order to assess the effect arising from Wilson loop autocorrelations, we 
form bins of the data of various sizes, though this reduces the amount of 
data available for estimating the covariance matrix even further. In fact, 
it turns out that the fits become unreliable beyond bin size 4 and before 
the binning error becomes stable. As a consequence we cannot exclude that 
the errors on ro/a are somewhat underestimated due to residual autocorrelations. 
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Mass dependence 

Our results for (r /a) are plotted in figure HI We note that the a\i q dependence 
appears to be rather weak, and hence we expect the data for the (purely gluonic) 
observable ro/a(afi q ) to be well described by polynomials of low order in afi q . In 
Table dU we collect the results obtained by fitting our data at different values of 
a\i q (see Table [TUT) to few simple functional forms, namely 

(I) : r /a + c 2 (afi q ) 2 , 
(II) : r /a + ci(a/i g ) , 
(III) : r /a + Ci(a// ? ) + c 2 (a/i 9 ) 2 . 

The ansatz (I) is inspired by the fact that with maximally twisted (unlike the 
case of untwisted) Wilson quarks the lattice fermionic determinant of the Nf = 2 
theory depends only quadratically on the bare quark mass. A weaker depen- 
dence on the bare quark mass can only appear via the effects of spontaneous 
chiral symmetry breaking on the static quark potential and would actually be a 
dependence Q on \afj, q \. This is the motivation for the fit ansatz (II), if it can be 
assumed that a/i g is sufficiently small to make the (a/i g ) 2 -dependence negligible, 
and (HI), if the (a/i q ) 2 -dependence is instead statistically significant. 

The fit based on the ansatz (I) describes our data rather well, as shown in 
figure HI suggesting that possible effects of spontaneous chiral symmetry breaking 
in the static potential at distances around 0.5 fm are negligible within our statis- 
tical errors. This interpretation is supported also by the other two fits: even if a 
/^-dependence of the type (II) cannot be ruled out completely, we observe that 
not only the xVd.o.f. of the fit (I) compared to (II) is better, but also the best-fit 
values of c 2 from fits (I) and (III) are more consistent between themselves (and 
less consistent with zero) than the best-fit values of c\ coming from fits (II) and 
(III). We would like to note that these findings are corroborated by analogous 
fits of the afi q dependence of the static potential at fixed values of r/a, i.e. in 
situations where no interpolation in r /a is involved. 



r /a 


ci x KT 2 


c 2 x 10" 4 


fit range 


X 2 /d.o.f. 


5.22(2) 




-0.08(2) 


B\ — B 5 


0.85 


5.22(3) 




-0.09(4) 


B\ — B A 


1.26 


5.28(3) 


-0.16(3) 




B\ — B 5 


1.10 


5.26(5) 


-0.12(6) 




B\ — Bi 


1.37 


5.22(8) 


-0.01(18) 


-0.08(9) 


B\ — B 5 


1.28 



Table 11: Results of the fits of the afj, q dependence of i'o/a according to the ansatz (I), (II) 
and (III) in the text. 



We conclude that the mass dependence is well described by the ansatz 
(I) and remark that an almost identical central value for r^/a at the chiral 

7 We are indebted to R. Sommer for very useful discussions on this point. 
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0.0001 0.0002 

(aji q ) 



Figure 4: Mass dependence of ro/a. The shaded area shows the error band of the quadratic fit 
(full line) to the data (circles). The additional plus symbols are further determinations of ro/a 
corresponding to different values of the fit parameters to the ansatz (|30[) . The spread provides 
an indication of the systematic error due to interpolation (see text) in ro/a. 



point is obtained from the ansatz (III), which also allows for a linear term 
in afi q . The ansatz (II) gives a central value for r Q /a at the chiral point 
lying two standard deviations above that from the ansatz (I). Finally we note 
that, if the ansatz (I) for the /independence of ro/a is used, the relative sta- 
tistical accuracy of our determination of r /a in the chiral limit is better than 1%. 



6 Some selected results 

In this section we present results for quantities related to the pseudoscalar (PS) 
channel. This includes, apart from charged and neutral PS masses and decay 
constants, also the renormalization constant Zy, which is specifically relevant to 
maximally twisted mass QCD. 
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6.1 Charged and neutral pseudoscalar masses 

Charged pseudoscalar meson mass 

To extract the charged PS mass mpg we consider the correlation functions 
discussed in sect. 12.11 We refer to this section and sect. 14.11 for a detailed discus- 
sion of how the correlation functions are evaluated and the errors are estimated. 
The results for the charged PS masses can be found in Table [121 

In order to make the effect of the longest runs at fi = 0.004 visible, we quote 
the results for run B\ a and the complete run B\ separately. While for S la we have 
1811 measurements made in the cyclic way explained in sect. 14.2} there are 895 
measurements for Bi performed moving the source time-slice randomly through 
the lattice. Even though in the latter case we have fewer measurements, they are 
more decorrelated because the single measurements are more separated in Monte 
Carlo time and because the distance of the position of the sources in Euclidean 
time for two consecutive measurements is on average larger. It is reassuring to 
see that results and errors are consistent between the two sets of data within 
errors. From this comparison it is also clear that moving the source time-slice 
randomly through the lattice is the most convenient of the two methods. 

In fig. [5] we show examples for effective masses in the PS channel at our 
lightest quark mass, extracted from the PS correlation function (with insertion 
of the d'-f^u operator) only. We plot the data for the three different choices of the 
interpolating operators, namely local-local, local-fuzzed, and fuzzed-fuzzed. One 
can see in fig. [5] that the three different operators give compatible results from 
t/a ~ 10 on. Hence we are confident that the ground state energy dominates for 
t/a > 9 and we chose the fit range accordingly. 

We also attempted to determine the energy of the first excited state of the PS 
meson from a 2-state fit to the 6x6 matrix of correlators. Even though we were 
unable to determine the first excited level in a reliable way from an unconstrained 
fit, fixing it to the theoretical value (3 times the ground state mass), as expected 
in the continuum limit, does allow an acceptable fit. 

This result is quite interesting, as with maximally twisted Wilson quarks one 
expects on general grounds also an 0(a 2 ) contamination from the 7r°(0)7r ± (0) two- 
pion state. Such a contamination becomes negligible, if compared to the expected 
three-pion state one in the continuum limit (taken at fixed quark mass). It should 
also be observed that when the pion mass is decreased, the two-pion contribution 
remains negligible with respect to the three-pion one until a 2 AQ CD <C e _mps *. 
For the range of mps and t values relevant for our data we find that two-pion 
contamination effects can hardly be detected despite our (small) statistical errors 
(see fig. [5]). 

Neutral pseudoscalar meson mass 

As discussed in sect. |3l the neutral PS meson can be created by interpo- 
lating fields that at maximal twist and in the twisted basis are of the form 
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Figure 5: Effective mass for the pscudoscalar channel from B\ lattice data. The effective 
masses obtained using 3 different interpolating operators as described in the text are shown. 



Run 


a/jL q 


amps 


afps 


ampcAc 






0.0040 


0.13587(68) 


0.06531(40) 


-0.00001(27) 


0.6114(85) 


B 1 


0.0040 


0.13623(65) 


0.06459(37) 


+0.00017(17) 


0.6136(19) 


B 2 


0.0064 


0.16937(36) 


0.07051(35) 


-0.00009(17) 


0.6096(21) 


B 3 


0.0085 


0.19403(50) 


0.07420(24) 


-0.00037(20) 


0.6115(22) 




0.0100 


0.21004(52) 


0.07591(40) 


-0.00097(26) 


0.6209(25) 


B 5 


0.0150 


0.25864(70) 


0.08307(34) 


-0.00145(42) 


0.6165(22) 



Table 12: Results for masses and decay constants in the charged pscudoscalar channel, PCAC 
quark mass and Zy . The results for the first three quantities come from a fit to a 4 x 4 submatrix 
with operators dj^u and idjoj^u, while for Zy we used the full 6x6 matrix. Note that the 
difference between the first two rows is just the length of the simulation. The time range of the 
fit was always 10 — 23 and the x 2 /d.o.f. was always smaller than one. 

XX an d X7o75X- We evaluate the correlator (both quark-connected and quark- 
disconnected pieces) with each of these operators at source and sink (also with 
local and fuzzed variants, thus giving a 4 x 4 matrix of correlators) as described 



U.J 



0.25 



0.2 



B °- 15 



0.1 



0.05 - 
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run 


a/jLq 


N 

1 > meas 


am PS,conn 


amp S 


a fps/ %a 


B 1 
B 3 


0.0040 
0.0085 


888 
249 


0.212(3) 

0.259(3) 


0.109(7) 
0.169(11) 


0.089(3) 
0.106(4) 



Table 13: Neutral pseudo scalar meson masses and decay constants at = 3.9 measured from 
every 10 trajectories at fi q = 0.004 and every 20 at \i q = 0.0085, as indicated; am,p S conn is the 
mass extracted from the quark-connected correlators only. 

above and in Appendix [B] We fit this correlator matrix to one or more states in 
the usual way. Based on our study of autocorrelations (see sect. (HI)), we compute 
statistical errors by a bootstrap analysis on blocked data where each block in- 
cludes measurements taken on configurations corresponding to a segment of MC 
history 80 trajectories long. 

Our results for the neutral PS meson are shown in Table [TBI Compared to 
ref. [TJ, we have increased statistics at fi q = 0.004 and we have employed the more 
refined fitting procedure explained above. In particular we used 4x4 matrix of 
correlators rather than a 2 x 2 matrix. We also include results at a second \i q 
value, fi q = 0.0085. In order to show the contribution of the quark-disconnected 
component to the neutral PS meson mass determination, we show appropriate 
ratios in fig. [61 

We have also evaluated the energies of neutral PS mesons with momentum 
2n/L (recall L/a = 24), obtaining (by use of the continuum dispersion relation 
E 2 = (27r/L) 2 + m 2 ) mass values consistent with those shown in Table [TBI 

The non-zero momentum results have the advantage that no vacuum subtrac- 
tion is needed for the neutral PS meson correlator and this provides a crosscheck 
of the approach we employed. For example, at fi q = 0.004 we obtain an energy 
of 0.309(27) which correspond to a mass 0.164^gQ in lattice units. 

Pseudoscalar meson mass splitting and related topics 

Concerning the PS meson mass, it is well known that with maximally twisted 
Wilson fermions, even in the theory with Nf = 2 degenerate quarks we consider 
here, there is difference of order a 2 (at fixed quark mass) between the neutral 
and the charged PS meson mass. Moreover the latter is very mildly affected by 
cutoff effects, once maximal twist is implemented in the optimal way of sect. 11.11 
as it follows from the formula mp s — mp s | cont = 0(a 2 fi q ) + 0(a 4 ) proved in 
refs. [TBI [T4"l 02]. Finally a lattice chiral perturbation theory analysis (see e.g. 
refs. |4*Bl [TBI |4*4"] ) shows that in the small fi q region the difference between the 
squared neutral and charged PS masses tends to an 0(a 2 ) quantity 

r 2 ((m° PS ) 2 -(m PS ) 2 )^c(a/r ) 2 , (31) 

which can be related to one coefficient (usually called C2) of the chiral effective 
Lagrangian of (twisted and untwisted) Wilson fermion lattice QCD. From our 
results we estimate c = —5.0(1.2) and c = —6.7(2.8) respectively at fi q = 0.004 
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Figure 6: Ratios of correlators of neutral to charged pscudoscalar (PS) mesons. The con- 
nected component of the neutral meson correlator is shown, as well as the total, which includes 
the disconnected contribution. The dotted curve illustrates the behaviour of the ratio of the 
correlator as obtained from our fitted values of the charged and neutral PS masses. 



and n q = 0.0085. These determinations are consistent within errors, as expected. 
Moreover, unlike in the case of lattice formulations with different gluonic ac- 
tions |45j . the sign of c is now in agreement with the first order phase transition 
scenario HH HZ] , as predicted by lattice chiral perturbation theory. 

Some results for the coupling of the neutral PS meson field to the divergence 
of the neutral axial current are also given in Table [13l The values of /p S are 
found to be quite close to those obtained for the charged PS meson (see fps in 
Table (EE])), if the estimate Z A = 0.76(2) pj EE] is employed. The differences 
between the central values are compatible with zero within 1 to 1.5 standard 
deviations of fp S (which is the least precise result of the two). 

This good consistency between the two channels is in striking contrast with 
the large cutoff effect present in m PS , which we observe as a large difference 
(about 50 MeV at a/iq = 0.0040 and a\i q = 0.0085) between its value and that of 
its charged counterpart, mpg. Theoretical arguments providing an explanation 
for the reasons why a large lattice artifact affects only m PS were presented in 
ref. (12] and will be further discussed in a forthcoming publication |50j . 
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6.2 PCAC mass 



As discussed in sect. (Il.ljl . we attempt to tune the value of mpcAc (see eq. (E])) 
to zero at our minimal fi q value, namely at afi q = 0.004. A definition of mpcAc 
equivalent to eq. ([3]) for time separations so large that the lowest PS meson state 
dominates, is given by 

m PS (0\A*\PS) , . 

mpcAC = —(o\p«\PS) ' a = 1 ' 2 - (32) 

These two matrix elements can be directly determined from a fit to the 4x4 
matrix involving the interpolating operators dj^u and id^^^u (or from the fit to 
the full 6x6 matrix, see sect. The results we obtain when a 4 x 4 matrix is 
used are summarised in Table [12] and shown as a horizontal line in fig. [7J It is 
important to notice that the condition discussed around eq. (flOl) is fulfilled for 
all our simulation points, and in particular for afi q = 0.004. 

In fig. [7J we also illustrate the time dependence of the local determination 
of the PCAC quark mass through eq. We see that the values of m PC Ac 

determined using eq. ([3]) in this way and eq. ( J32l) agree very well between them- 
selves, in the t-region where the ground state pseudoscalar meson dominates, as 
expected. 

Compared to ref. [JJ we now have a result for ampcAc available for the large 
statistics run B\. It is reassuring that there is full consistency between the 5000 
trajectory run and the run extended up to 10000 trajectories. This makes us 
confident that our error estimate is realistic. Our results for ampcAc & s a function 
of the bare quark mass afi q are illustrated in fig. El where results from the full 
ensemble, B%, and those from the smaller ensemble, Bi a , are separately shown. 

6.3 Pseudoscalar decay constant and Zy 

Using the exact lattice (twisted basis) PCVC relational 

(d;v;(x)O(0)) = -2^ ab (P b (x)O(0)) a = 1,2, (33) 

where <9* is the lattice backward derivative, O a local lattice operator and V^{x) 
the 1-point-split vector current 



1 

4 

X{x + pL)T a Ul{x){ lil + r)x»{x)] , (34) 



= ~ A [xi^TaU^x)^^ - r) X {x + A) + 



we can also compute (in the charged meson channel) the pseudoscalar meson 
decay constant with no need of any renormalization constant (see [2} [HUH]) from 



We recall that at maximal twist the twisted basis vector current corresponds to the 

lb 



axial current e 3ba A'^ (a,b = 1,2) in the physical quark basis 
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Figure 7: PCAC quark masses versus time separation (from eq.0) for three values of the mass 
parameter, \i q = 0.004, fi q = 0.0064 and fj, q = 0.0085. The lines represent the central values of 
the PCAC quark mass from a fit to a 4 x 4 matrix of observables for a i-range 10-23. 



the formula 



/] 



PS 



2/i q 



m 



\(0\P a \PS)\ 



1,2. 



PS 



Based on the exact relation (1331) and noting that 



0(a 2 



(35) 



(36) 



where V® is the naive current defined in eq. (1471) and <9 M the symmetric lattice 
derivative, we can determine the (scale-independent) renormalization constant 
Zy through (a, b = 1, 2 and Xq ^ 0) 



Z v (a/j, q 



2/^Ex£ 3 ( p (^) P (°)) 
E x 9 (^(x)P fc (0)) 



and 



lim Zy(a/jb q ) 



(37) 



with only 0(a 2 ) cutoff effects. The results for a/ps extracted using eq. ( !35l) are 
collected in Table [12 and will be discussed further in the next section. 

In the same Table we quote the results for the determination of Zy that 
is based on fitting the large time behaviour of eq. (1371) where the pion state 
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Figure 8: The PCAC quark mass ampcAC as a function of a/i g 
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In practice the relevant matrix elements are extracted from a factorising fit (see 
eq. ( fl4l) ) to a (4 x 4) matrix of correlators, with entries given by the expectation 
values of £ x P 1 (x)P 1 (0) and £ x V 2 (x)P 1 (0) with local-local and local-fuzzed 
operators. 

The an q dependence of Zy is everywhere very weak (see fig. [H]), which makes 
the extrapolation to the chiral point easy and almost irrelevant, and leads to a 
rather precise result: Zy = 0.615(5). The quoted error is a conservative estimate 
of the total uncertainty on Zy, inferred from the data in the last column of 
Table [T2] and their statistical errors. 

Another determination of Zy(an q ) is obtained by evaluating the ratios of cor- 
relators in the r.h.s. of eq. (37) for a number of time separations (xo/a > 10) and 
taking the average over them. This alternative method to evaluate Zy has been 
presented in ref. jlH]- Applying it, for instance, to the available (~ 900) correla- 
tors computed on the whole ensemble P>\ yields the result Zy(afi q = 0.0040) = 
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0.62 




0.61 - 



Figure 9: Extrapolation of Zy to zero quark mass at (3 = 3.9. The data are consistent within 
errors with the constant behaviour shown in the figure. 



0.6101(2). Both approaches, the one discussed in some detail above and the one 
of ref . [IE] , provide precise determinations for Z v and the virtues of both methods 
will be further discussed in ref. |49j . 

7 Chiral Perturbation Theory analysis of /ps 

and mps 

In this section we present the details of the comparison of our data with Chiral 
Perturbation Theory (xPT) predictions. The main results have been already pub- 
lished in ref. pp. Here we provide further information about our fitting procedure 
and error analysis. The main goal of this section is to explain the fitting proce- 
dure and error determination using chiral perturbation theory. We will therefore 
only use the ensembles B\-Br )1 i.e. the ensembles that have already been discussed 
in ref. p]. For this limited set of data at only one value of the lattice spacing of 
a = 0.087fm, we are not sensitive to higher order effects of chiral perturbation 
theory. We will therefore restrict ourselves to a 1-loop analysis of the data on M PS 
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and fps- Nevertheless, we use the 2-loop chiral perturbation theory expressions 
and vary parameters of the corresponding formulae to see the possible effects, 
if the 2-loop order would be important to describe our data. In particular, this 
2-loop investigation confirms that the here chosen dataset is indeed not sensitive 
to higher loop corrections. 

Some preliminary results at larger volumes (L/a = 32) and finer lattice spac- 
ing (j3 = 4.05) have been already presented [TJIQI)]. However, the present work is 
focused on the details of the analysis of the data points presented in ref. pQ. The 
study of the volume and scaling dependence will be presented elsewhere. Notice, 
however, that, w.r.t. ref. pQ we have a larger statistics at the smallest quark mass. 

Our raw data for ampg and a/ps are determined as described in sect. 14.11 
Results are reported in Table [T2J As said above, there is no need to compute 
any renormalisation constant in order to make contact with the corresponding 
physical quantities. 

In our ;\TT based analysis we have to take into account finite size corrections 
because on our lattices at the lowest and next-to-lowest /! g -values they turn out 
to affect amps and, especially, a/ps in a significant way. 

We have used continuum xPT to describe consistently the dependence of the 
data both on the finite spatial size (L) and on the bare quark mass This 
might look inappropriate in view of the existence of a large additive 0(a 2 ) artifact 
in the neutral pion mass squared 0. However this is not so, because theoretical 
analyses carried out in the framework of lattice xPT [T3] and Symanzik expan- 
sion (complemented with soft pion theorems in the continuum theory) WI\ 
show that, if maximal twist is implemented as discussed in sect. 17.21 the charged 
pion squared mass differs from its continuum counterpart only by 0(a 2 /i) and 
0(a 4 ) terms, while the charged pion decay constant is affected by (chirally non- 
enhanced) discretization errors of order a 2 . Moreover the Symanzik expansion 
analysis is applicable for all spatial volumes L 3 , provided L is large enough to 
justify the use of soft pion theorems in the continuum theory at the quark mass 
of interest. This entails that also the L-dependence of the charged pion squared 
mass and decay constant is expected to be essentially continuum-like. These ex- 
pectations are also supported by preliminary and still partial results we obtain 
at different lattice resolutions and different physical volumes [T5"| [T5] . Last but 
not least, the continuum xPT formulae appear to describe well our data, as we 
are going to show below. 

We fit the appropriate (Nf = 2) next-to-leading-order (NLO) xPT formu- 

9 Theoretical arguments have been presented [33] suggesting that this lattice artifact is an 
exceptional, though important, case, because it is related to the large value of a continuum ma- 
trix element appearing in the Symanzik expansion of the 7r°-mass and does not stem from large 
coefficients multiplying the dimension five and six terms of the Symanzik effective Lagrangian. 
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lae [521 G3] 



m ls( L ) 



2B fi q 



[l+flog(25 /z,/A 



(39) 



fp S (L) = F [1 — 2^i (A)] [1 - 2eiog(2B ^/At)] , 
to our raw data for mps and fps simultaneously. Here 

f = 2 J B /i g /(4vrF) 2 , A = y/2B of i q L* . 



(40) 



(41) 



The finite size correction function gi(X) was first computed by Gasser and 
Leutwyler in ref. [52] and is also discussed in ref. [53] from which we take our 
notation (except that our normalisation of f n is 130.7 MeV). In eqs. (I3"9"]) and (1401) 
next-to-next-to-leading order (NNLO) corrections are assumed to be neg- 
ligible (this assumption is critically discussed in sect. 17.11) . The formulae above 
depend on four unknown parameters, Bq, F, A3 and A4, which will be determined 
by fitting to our data. 



7.1 Statistical errors 

In order to estimate the errors on the fit parameters it is important to account 
for both autocorrelation and cross-correlation of the data. We have exploited two 
different methods to do so. 

Method A 

The first method (see also [TU] ) consists in computing the full covariance matrix 
of our data for afps and (ampg) 2 and include it in the computation of \ 2 

^ = Y f {\H-Fi)Vry{y j -F j ), (42) 

where V is the covariance matrix 

V id = cov(j/i, yj) = cov((^ - F i ){y j - Fj)) , (43) 

normalised so that the diagonal elements coincide with the squared standard 
error, and we have expressed the ansatz, eqs. ( 139]) and ( 140]) . in the form 

y i = F i (x,6). (44) 

10 We stress that £ defined here should not be confused with the continuum matrix element 
£,r introduced in eq. (|5]>. 

11 As we have data from independent simulations (ensembles B\ to -B5), in the present case 
the covariance matrix will be block-diagonal with five blocks. 
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Here we denote by y^ the primary measured quantities (in this section: yi = 
(am PS \ Bl ) 2 , y 2 = af PS \ Bl , y 3 = (am PS \ B2 ) 2 , ...), by x the independent (error- 
free) variables (in this section: only ) and by 6 the parameters to be 
determined by the fit (here: 9\ = 2aBo, 62 = aF, #34 = log(a 2 A| 4 )). The error 
on the parameters are thus given by 

(A0 a f = (Vg^V-'Ve^)- 1 . (45) 

The autocorrelations of (ampg) 2 and af P $ have been estimated both by data- 
blocking and by means of the T-method, as discussed in sect. 14.11 Both ap- 
proaches indicate (see Table [3]) that by combining data into blocks of 32 mea- 
surements each (this corresponds always to more than 60 MC trajectories) the 
resulting blocked data are safely uncorrelated. These blocked data are thus used 
to evaluate the covariance matrix, the y 2 and the errors on the fit parameters as 
discussed above (eqs. H2HH]). In this way the possible effect of cross-correlations 
among the observables is included in the covariance matrix and therefore properly 
accounted for in the fit procedure. 

In some of the checks that we are going to present below it will not always 
be possible to reduce the xPT formulae to the form of eq. (1441) . This happens 
in particular in the following cases: when computing directly a/ps as a func- 
tion of amps, when including the effects due to a non- vanishing ampcAc, when 
including higher orders in Finite Size Effects (FSE) calculations, as computed 
in [53], or when we will eventually study the scaling dependence for different 
lattice spacings a. In all these cases the xPT formulae can be expressed in the 
more general form Gi(y,x,6) = and the errors are given by the formula |54j : 
(A6 a ) 2 = iye a G T {y y GVW y G T )' l We a G)- 1 . These errors are obtained as prop- 
agation from the known errors on x, therefore they do not depend on how good 
the fit is. The quality of the fit is expressed as usual by the quantity x 2 /d.o.f.. 

To provide a further check of possible effects of cross-correlation, the fits were 
also performed by dividing the data set into two subgroups each of half the size. 
The data for ampg were taken from one subgroup of gauge configurations and 
those for af P s from the other, ensuring in this way absence of cross-correlation. 
Errors scale as \/2, i.e. as expected from halving the statistic, which indicates 
a negligible effect of cross-correlations in the full data set. Stability was also 
checked against different choices of subgroups. This result is confirmed by the 
observation that if we suppress the off-diagonal terms in the covariance matrix, 
our error bars are affected only at the percent level. 

Method B 

Method A is standard, and of course unbiased if for all observables the data are 
distributed in a Gaussian way (which we checked explicitly to be the case to 
a good approximation) and if the functions F (or G) have a sufficiently linear 
behaviour around the relevant values of their arguments. An even safer estimate 
of the final errors can be obtained with the bootstrap method [3B1 ES] • 
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Method A Method B 



2aB Q 5.04(7) 5.04(7) 

aF 0.0522(7) 0.0522(7) 

log(a 2 A 2 ) -1.90(11) -1.91(10) 

log(a 2 A 2 ) -1.00(4) -1.00(4) 

X 2 /d.o.f. 1.0/4 0.9/4 

Q 0.91 0.92 



Table 14: Comparison of fit results from methods A and B. 

To apply the bootstrap analysis method to our data set we proceed as follows. 
In order to account for autocorrelations we first form bins of 32 gauge configura- 
tions for each value of \i q . Out of the blocked data we generate 1000 bootstrap 
samples. The size of each sample is chosen as large as the full (blocked) data set. 
From the 1000 bootstrap samples we obtained 1000 observations for 2aB , aF, 
log(a 2 A 2 4 ), and " 3i 4 = log(A 2 i4 /m 2 ), respectively. Error estimates are then com- 
puted as prescribed by the bootstrap method, i.e. by the standard deviation over 
the (equally weighted) 1000 samples. Incidentally we remark that this procedure 
takes the cross correlation between amp s and a/ps correctly into account. In the 
1000 fits we performed we have hence always used only the diagonal elements 
of the covariance matrix (fixed to their central values, i.e. to the square of the 
statistical errors on amp s and a/ps, see Table [12]) as weights to evaluate the \ 2 
formula ( 1421) . In this specific application of the bootstrap method the errorbars 
on the basic quantities amp s and a/ps are still needed, since the observables of 
interest, the low energy constants (LEC), are defined through minimization of 
the x 2 °f the simultaneous fit to the formulae eqs. (131^ - (14*01 . Note that to 
safely employ the bootstrap method data need not have a Gaussian distribution, 
nor do the constraints, defined by the formulae, need to be linear. The 
bootstrap method may become expensive if single fits are significantly computer 
time demanding. 

Both methods A and B give consistent results, as shown in Table [TH In this 
paper we use the same setup as in ref. p], but we employ a somewhat larger 
statistics. The results are consistent. In addition to the error estimates we quote 
the value of x 2 an d the merit figure of the fit defined via 

Q = 1 - P( X 2 /2, d.o.f./2) , 

where P is the incomplete Gamma function |56j . 

7.2 Discussion of systematic errors 

The error bars quoted in Table [TH are only statistical. As we also stressed in 
ref. [TJ, a number of systematic effects are expected. Here we present some 
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including q/i 9 = 0.015 m q = ^ (Z A m PC ac) 2 + /J% 

2oBq 5.06(5) 5.05(6) 

aF 0.0508(5) 0.0521(7) 

log(a 2 A|) -1.93(6) -1.87(11) 

log(a 2 A 2 ) -0.89(2) -0.99(4) 

X 2 /d.o.f. 10.3/6 0.55/4 

Q 0.11 0.97 



Table 15: Comparison of fit results from different setups, as explained in the text. 



checks we performed in order to estimate the actual magnitude of these systematic 
effects. 

As a first, simple check on the possible impact of neglected NNLO terms 
on the results presented in Table HH we have also included the heaviest point 
(the one at a\x q = 0.0150) in the standard fit to the formulae fl3U|) - fj4*0l . In this 
case the results are still compatible with those in Table [14] within 1.7 standard 
deviations, but the % 2 /d.o.f. of the fit jumps from 0.24 to 1.7. This increase of 
X 2 /d.o.f. is mainly due to the point at a/ps at a\i q = 0.0150, as we noted already 
in ref. pQ. The results of the fit are displayed in the second column of Table [Tol 
This suggests that only the first four quark mass points should be used when 
comparing our data for a/ps and am-ps with NLO xPT, as was done in ref. pQ. 

It is also very interesting to see how much the tiny deviations from maximal 
twist corresponding to the (statistically compatible with zero) measured central 
values of mpcAc affect our results for the low energy constants discussed in this 
section. To address this question we introduce the definition of bare quark mass, 

m q = ^/(Z^mpcAc) 2 + Z^ 2 ) which holds for generic twist angle up to neglected 

O(ampcAc) and 0(a 2 ) terms. Moreover, in order to take into account the axial-T3 
transformation properties of the current entering the formal definition of /ps, at 
the same level of accuracy, the value of a/ps should be corrected into af-p$m q / n q . 
We remark that this is obtained automatically if fps is evaluated from eq. ( 1351) 
with fi q replaced by m q - this can be related to the invariance of the operator 
P 1 ' 2 , a matrix element of which appears on the r.h.s. of eq. (|3"oT) . under axial- t% 
rotations. The results of this analysis, where we set Za = 0.76(2), as found at 
(3 = 3.9 in ref. [H], are shown in the last column of Table fT5l It is reassuring to 
see that, thanks to the good precision we could reach in setting mpcAc to zero, 
the low energy constants of interest here are left essentially unaffected by this 
kind of correction. 

We now consider the finite size corrections. In ref. [T] we estimated them 
with the help of the formulae of ref. [52]. A nice feature of these formulae is that 
they introduce no new parameter. However, they are only the first term of an 
expansion. Hence, the question is: how large is the residual uncertainty in FSE 
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a\i q 


lo [52J 


lo [S 


a 


nlo [1 


31 


nnlo [53] 


0.0040 


0.64 % 


0.42 


% 


0.50 


% 


0.21 % 


0.0064 


0.29 % 


0.16 


% 


0.21 


% 


0.10 % 


0.0085 


0.16 % 


0.08 


% 


0.12 


% 


0.06 % 


0.0100 


0.11 % 


0.05 


% 


0.08 


% 


0.04 % 



Table 16: Percent Finite Size deviation (mps(i) — mps(oo))/mps(oo) predicted by xPT for 
our data points. Note that nlo and nnlo include only the last order and not the previous one(s). 
According to ref . |57j , comparing nlo and nnlo (not lo and nlo) gives a reliable indication about 
the convergence of the expansion. 

due to this truncation? To go beyond the first term in the framework of ref. [52] 
is difficult. For the pseudo scalar mass the FSE corrections at two loops in xPT 
have been computed in ref. [S3]- However, one can do better using the kind of 
%PT expansion suggested in ref. [57], for which results are also given in ref. [53J. 

With the help of the results from ref. [53] we can assess the stability of the 
prediction both by comparing the two approaches and by studying the conver- 
gence of the expansion of refs. [571 [53] . One should also notice that higher orders 
do introduce new parameters. Since it is not realistic to fit them, we will instead 
look at the stability of the prediction while changing those parameters in a "rea- 
sonable" range. The "reasonable" range is suggested in ref. [53] and is based on 
phenomenological grounds. 

To avoid confusion, we remark that the results of ref. [52] are given as an 
expansion in powers of 1/Fq, while ref. [53] uses an expansion in 1/F n . This is 
the only reason why the first term of ref. [53] does not coincide with ref. [52] • 

In Tables [TBI and [T7I we show the percent deviation obtained using the formulae 
from refs. [52] and [53] at different orders. Note that the new low energy constants 
(LECs) that at higher orders of xPT are relevant for FSE are fixed to their central 
values estimated in ref. [53]. See the comment below about their impact. To 
distinguish the expansion of the FSE effects from the usual xPT expansion we 
will use a lower case notation (lo, nlo, nnlo) to denote the former one. The two 
expansions are of course related, but since the FSE also depend on the lattice 
size L, there is no reason to truncate the chiral expansion for FSE at the same 
order as the usual %PT expansion. Here, for instance, we will use the NLO %PT 
formulae, but we will compare FSE at lo, nlo and nnlo. 

The convergence of the FSE expansion is expected to be good for all our 
data points since the smallest value of mpsL is larger than 3. We recall that, 
according to ref. [57], the comparison of lo and nlo is not a good indicator of 
the convergence of the expansion. This should be rather checked by comparing 
nlo and nnlo. According to all our estimates only the FSE at the lightest point 
{afjiq = 0.004) are relevant, while those at larger quark masses are always smaller 
than statistical errors. For instance, the deviations in mps are barely larger 
than its statistical errors (which amount to about 0.5%). In order to check the 
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afi q 


lo [52J 


lo [53j 


nlo [53J 


0.0040 


-2.57 % 


-1.68 % 


-0.76 % 


0.0064 


-1.15 % 


-0.63 % 


-0.30 % 


0.0085 


-0.64 % 


-0.32 % 


-0.16 % 


0.0100 


-0.44 % 


-0.21 % 


-0.11 % 



Table 17: Same as in TableQH but for (f PS (L) - /ps(oo))// ps (oo). 



dependence of the predicted FSE corrections on the LECs entering only at nlo, we 
changed randomly the value of the latter within the "reasonable" range suggested 
in ref. [53]. We saw that nlo and nnlo FSE corrections are affected only at the 
level of about 20% (lo corrections are obviously unaffected) by such changes. 

Up to this point we have only considered the %PT at NLO (however correc- 
tions as high as nnlo are included in FSE calculations) implicitly assuming that 
NNLO contributions are negligible. This is reasonable, since formulae with 
only NLO corrections yield a very good fit of the data at the lightest four quark 
masses, in spite of the fact that the expansion parameter, £ = 2Bofi q /(4:7rF) 2 , 
is not always very small. It is thus important to assess how much NNLO terms 
would affect our results. 

The NNLO corrections relevant for mpg and /pg have been calculated in 
ref. [58]. Here we use an expression which is easier to compare with lattice data, 
namely the one of refs. [EJ [60] which reads 
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where £ = 2B fi q /(4nF) 2 as before, M 2 
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It is not realistic to attempt a fit of all the coefficients involved in the full 
NNLO expressions at least with the limited set of data used here. Rather we 
fix the parameters Ai, A 2 , kp and ku to the values suggested in ref. [53]. Since 
no estimate for ku,F is available, we take ku,F = 0. Redoing the fit in these 
conditions we can check how much NNLO terms change the results of Table [141 
The new fit results are shown in the second column of Table [TH1 In order to 
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NNLO as in (£3 


5Ai = ±33% 


5A 2 = ±5% 


k M = ±1 


kp = ±1 


2aB 


4.80(6) 


-0.66% 


-0.20% 


3.2% 


0.07% 






3.44% 


0.26% 


-2.5% 


-0.12% 


aF 


0.0536(6) 


0.60% 


0.16% 


-0.19% 


1.9% 






-1.7% 


-0.19% 


0.21% 


-2.1% 


log(a 2 A|) 


-2.13(12) 


-9.6% 


-1.2% 


-29% 


-1.3% 






-5.9% 


0.87% 


26% 


1.5% 


log(a 2 A 2 ) 


-1.00(5) 


-4.6% 


-0.50% 


1.3% 


24% 






-0.35% 


0.34% 


-1.3% 


-26% 


xVd.o.f. 


0.085 


1.7 


1.1 


0.48 


1.4 






0.15 


0.82 


1.8 


0.73 



Table 18: Fit results, including NNLO xPT. The second column shows the results obtained 
with the choice of A12 suggested in [53] and kM,F = 0. The other columns give the percent 
correction due to changing the corresponding parameter in the indicated range. For each line, 
the upper (lower) number corresponds to the higher (smaller) boundary value of the interval. 



further estimate to which extent these numbers are sensitive to a change in the 
parameters which were held fixed, we decided to change them one by one within 
the range proposed in refs. EZj, and perform a new fit for each one of these 
values. As for ku and kp, it is difficult to tell what is a reasonable range, since, 
as we said, no estimate is available for them. On general grounds the values of 
ku,F are expected to be of 0(1) and somewhat arbitrarily we assume a variability 
range ku,F = ±1- This choice is also justified by the fact that larger variations 
quickly lead to very bad x 2 - The results of this elaborated procedure are shown 
in columns 3 to 6 of Table [TBI Most effects are not significant if compared to 
statistical errors, as they are never larger than a few standard deviations. It 
should be noted, however, that A3 appears to be rather sensitive to kM and 
similarly A 4 to kp. These LECs can deviate by about 25% when setting ku,F to 
+ 1 or to —1. We mention that changes of the LECs similar to those reported in 
Table HB] are also obtained if the NNLO terms in eq. (j4~6l) are replaced by simple 
polynomial terms, like Pm,f£, 2 (with no logarithms), and the free parameters Pm,f 
are set to their best fit values. 

7.3 Comments 

In summary the discussion developed in this section shows that at least the 
systematic errors coming from the unknown NNLO terms involving /cm,f may be 
significantly larger than the statistical ones, mostly because the adopted range of 
values was, to some degree, arbitrarily chosen. However, as already said above, 
using only the datasets Bi-B 5 a reliable estimate of systematic uncertainties on 
B, F, A 3 and A 4 from the NNLO corrections is not possible. A better assessment 
about the magnitude of NNLO effects will be attempted elsewhere [HI] using 
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ETMC data at different lattice spacings. 

Although FSE to our simulation data turn out to be less than a few percent, 
we have made a special effort to compute them quite accurately, because their 
impact on LECs cannot be neglected, as their magnitude is comparable to the 
size of our statistical errors. The computation of FSE made in ref. [53] represents 
a considerable improvement on the classical estimate of ref. [S2| , as uncertainties 
on the extra LECs entering the former computation at high orders have little 
impact on the results. Actually, the validity of the predictions of FSE from %PT 
can be checked by performing simulations on lattices of increasing size in physical 
units. Preliminary results have been presented in ref. [T5] . 

8 Summary 

In this paper we have illustrated and discussed a number of details concerning 
unquenched simulations of Nf = 2 mass degenerate Wilson quarks at maximal 
twist. We have explained in sect. [Hour criterion on how to tune the theory to 
maximal twist. In particular, we provided theoretical arguments for our choice of 
m PCAc/V<? < 0.1 and showed that an error AmpcAc/V<? < 0.1 is appropriate for 
this purpose. Useful formulae for quark bilinears and their physical interpretation 
in different quark bases (twisted and physical) are collected in Appendix [A] 

We have then discussed in sect. [2] the methods we have used to compute 
charged meson correlators emphasizing the effectiveness of employing (fuzzed) 
stochastic time-slice sources in the so-called "one-end trick". We have demon- 
strated that this method complemented by a random choice of the source location 
leads to a significant noise reduction, at least for two-point correlators in the me- 
son sector. 

The computation of neutral mesons and, in particular, quark-disconnected 
contributions has been described in sect. [3] and in the corresponding Appendix [Bi 
We have spelled out the reasons for using stochastic volume sources which can 
be employed in combination with efficient variance reduction methods. All these 
technical improvements have allowed us to compute quark-disconnected contribu- 
tions on our sets of unquenched gauge configurations to an acceptable accuracy. 

In sect. H]we have illustrated the main features of the MC algorithms used in 
our simulations showing that the resulting autocorrelation times are small enough 
to allow for a trustworthy error analysis of physical observables. We also explain 
how our error analysis of the data was performed owing to the use of T- and 
binning- methods. 

The force parameter r can serve as an important physical quantity to check 
the scaling behaviour towards the continuum limit. We have provided in sect. [5] 
a comprehensive discussion of the methods we have used to extract tq on our 
configurations. It turns out that with the present data an accuracy of better 
than 1% can be reached for r in the chiral limit. It is also found that r has a 
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mild quark mass dependence which is consistent with being quadratic in fi q . 

Various results for the charged and neutral pseudoscalar masses, the untwisted 
PCAC quark mass and the renormalization constant Zy are collected in sect. [6j 
In particular, we show "effective mass" plots demonstrating the stability of the 
Euclidean time plateaux, which enables us to extract precise results for mesonic 
quantities. 

Finally, we have detailed in sect. [7] how our analysis of the data on 

mps and fp$ has been carried out, explaining how we get errors on the fitted 
low energy constants of the effective chiral Lagrangian, Bo, F, A3 and A4. In 
addition, we have analyzed the effects of higher orders in x?T on the stability of 
fit parameters and discussed the finite size effects. 

We consider the present paper as a technical reference work of our collabo- 
ration. The methods described here have been and will be extensively used in 
our ongoing future research on lattice QCD employing maximally twisted Wilson 
fermions. 
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Appendices 



A Quark bilinear operators in the twisted basis 

We give in this appendix the expression of a number of bare quark bilinear oper- 
ators that are relevant for the topics of this paper. The operators are expressed 
in terms of i) simple composite fields (recall 75 = 70717273 and = i/2[y fl , 7„]) 
in the twisted quark basis, where the fermionic action takes the form (jTJ) , 

S°(x) = x(x)x(x), P a (x) = xOhsyXO), 

A fa) = x(z)7M75yX(z), V£ (x) = x{x)l^x{x), (47) 

T »v( x ) = X(x)a^—x(x), T° u (x) = x(x)<r^x(x) . 

and ii) the twist angle u, where tana; = fi q /(m — m crit ) and am cr i t is determined 
as discussed in section 11.11 The expressions we get are 

jcos(u)A- + e^sm(u;)V^ (a = 1,2), 
M "U (« = 3), (48) 

JcosM^ + e 3 ^sin(u;)^ (a = 1,2), 
V »~U (« = 3), ( ) 



(50) 



cos(cu)P 3 + i|sin(cj) 1 S° (a = 3), 
P a (a = 1,2), 

S'° = cos(u)S° + 2i sm{u)P 3 , (51) 



1 (a 1, 2), 



cos 



HT 3 ,-zie^sinHT A ° p (a = 3) 



These expressions follow from the relation between twisted basis (x) and physical 
basis (i/j) quark fields, which (see eq.(pQ) and ref. [3]) reads 

X = e-^' 2 ^ , x = ^e~ 1 ^ 12 , (53) 

and the (obvious) definitions of the bare primed operators in terms of physical 
basis quark fields (a = 1,2,3) 

A'; = <K*) 7M 75y^), C = Hxh^(x), 

P' a = ^(x) 7 5y^), T£{x) = ^(x)a wt y^(x), ( 54 ) 
/o 



S m = ip(x)tp(x) . 
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All these bare operators renormalize multiplicatively, with the exceptions of P' 3 
and 5"°, which undergo an additive mixing with the identity (cubically divergent 
for P' 3 , quadratically divergent and vanishing as fi q — > for S'°). For the ex- 
pression of renormalization constants as functions of uj and the renormalization 
constants of standard Wilson quark bilinears and further details, see ref. [3]. It 
should be remarked that substantial simplifications occur for u = ±7r/2 (max- 
imal twist) in almost all formulae above. Moreover at maximal twist also the 
formulae for renormalization constants [3] get much simpler than at generic u. 

B Evaluation of disconnected loops 

The quark-disconnected (simply "disconnected" in the following for brevity) com- 
ponents of correlators are intrinsically noisier than the connected components, so 
it is essential to evaluate them as accurately as possible. For this purpose we 
need to compute the disconnected loops at every t value and for as many gauge 
configurations as possible. This can be achieved by using the stochastic source 
methods as we now discuss. The goal of the approach is to have an error aris- 
ing from the stochastic nature of the method which is smaller than the intrinsic 
variability associated with varying t and gauge configuration. If this is achieved, 
then the stochastic error is negligible in the sense that any further improvement 
in the signal can only be obtained if more gauge configurations are employed. 

As discussed in sect. 12.11 the basic idea is to use stochastic sources (£ ) having in 
general support on the whole lattice and solve the linear system for the quantities 

(f> = M- 1 £, (55) 

where M is the lattice Dirac matrix for a given flavour. The equation above is 
the same as eq. ( fl6l) . with the omission of the noise sample label r (to lighten 
notation). Note also that in this appendix the normalization of M is taken such 
that, if -Di att denotes the two-flavour Dirac matrix in eq. ((T]), then 

M u = 2«tr[aA«*t(l + t 3 )/2] = A + H, A = l + 2k W75 , (56) 
with H the usual Wilson first-neighbour hopping matrix. It follows that 

Y^X<f>\ R = J2 XM ~ X + noise (57) 

where the symbol [■■■]r refers (as in sect. 12. ip to the average over R samples 
of the stochastic source, the symbol denotes the sum over colour, spin and 
space-time indices and X can be (almost) any structure we wish to evaluate, like 
7- matrix, gauge links, Fourier factor, cos(kx), etc... It should be observed that in 
evaluating the disconnected contributions to the neutral meson correlators each 
one of the two quark loops arising from Wick contractions must be averaged over 
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completely independent samples of stochastic sources for the purpose of avoiding 
unwanted biases. Moreover, for each quark loop diagram, the sum in eq. (1571) 
is restricted to one single time-slice, while still ranging over all color, spin and 
space indices. 

A method we employed to reduce the variance of the stochastic noise with- 
out much additional computational effort is the hopping-parameter method [25J. 
This relies on the observations that the first four terms in the hopping parameter 
expansion of ^2 XM~ X can be easily evaluated exactly on each gauge configu- 
ration and that replacing their stochastic estimates with the exact values sig- 
nificantly reduces the variance. In fact, writing M u (see eq. ( 1561) ) in the form 
M u = (1 + HB)A, where B = 1/A, one easily obtains the identity 

1/M U = B- BHB + B(HB) 2 - B(HBf + {l/M u )(HBf , (58) 

which can be used to give 

^X/M u = ^{X(B - BHB + B{HB) 2 - B{HBf + (1/M U )(HB) 4 )} . 

(59) 

The last term in eq. (I5TJ1) can be evaluated stochastically because 

X(1/M U )(HB) 4 = lim [e{HB) A X^)} R (60) 

it— >oo 

Since = 75-^75 and 75 commutes with B, the last formula can be rewritten 
in the form 

X(1/M U )(HB) 4 = \im[( l5 (B^H)%0*X(f)} R . (61) 

Thus four extra multiplications of the source £ by B^H are needed. This is 
a negligible overhead compared to the inversion needed to obtain <fi. The first 
four terms in eq. (I5D1) do not involve 1/M U and can be, as said above, evaluated 
straightforwardly for any choice of X. For a local operator X, the only non-zero 
contributions are from the first term if X is proportional to 1 or 75 and the third 
term if X is proportional to 75. For a non-local operator X whose length of 
spatial path is more than four lattice hops (as used in this paper), the first four 
terms are all zero. 

This variance reduction method reduces the standard error of the stochastic 
samples by a factor of 1.5 or more in our case. This is valuable (it saves a factor 
2-3 in computational time), but for twisted mass QCD a much more powerful 
method is also available, although it applies only to the case J2X(1/M U — 1/Md). 
This last method can be, and has been indeed, used in many important applica- 
tions, essentially all those where one has to evaluate correlators with insertions 
of neutral meson operators of the form (in the twisted basis) x^ T 3Xi with any 
Dirac matrix T and the flavour matrix T3. Interpolating fields of this type occur 
e.g. in the two-point correlators for 7/, fo and tt° mesons (actually only one of the 
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possible operators for 7r°), as one can see from the table for the (twisted basis) 
neutral meson operators reported in sect. [31 

The powerful method alluded above relies on combining the identities 

(M d -M u ) = -4ma/i 9 7 5 (62) 

and 

(l/M d )(M d - M u )(l/M u ) = 1/M U - 1/M d (63) 

to get 

1/M U - l/M d = -AiKa N (l/M d ) l5 (l/M u ) . (64) 

The latter relation already serves as a method of variance reduction because 
the explicit (small, in our simulations) factor of a\x q reduces the magnitude of 
the fluctuations. On top of that, an even more important point is that the 
r.h.s. of eq. fl64l) can be evaluated very effectively with the help of the "one- 
end-trick" [T91 |2"U] and no further inversions. In fact, since M\ = j^M^^, one 
has 

J2X(1/M U - l/M d ) = -AiKa^Y.X^il/M^il/Mu) , (65) 
which can be evaluated with noise/signal ratio of 0(1) via 

X{l/M u - l/M d ) = -AiKafi g ^[0*X 75 0] R + noise , (66) 

where (we recall) <p — 0-/M u )£ and <p* = (1 / 'M u y . Apart from the explicit sum 
denoted by J^, the r.h.s. of this formula contains an implicit sum over the space- 
time indices of the stochastic source £ in and 0*, which contributes to reduce the 
variance as it receives contributions from the whole lattice (space-time) volume. 

To give an idea of the effectiveness of the method based on eq. (I6"6"j) we consider, 
as an example, the special case X = where one obtains 

^ 75 (1/M U - 1/M d ) = 4KCLfi q Y^W^R + noise . (67) 

At P = 3.9 and fi q = 0.004 (ensemble B±) the method based on eq. ( 1671) yields 
an error which turns out to be 6 times smaller than what would be obtained 
with a conventional stochastic volume source. From the measured stochastic 
contribution to the signal, as well as the observed total fluctuation, one can 
extract the intrinsic variation stemming from the statistical fluctuations of the 
gauge field. The goal of the stochastic method is to have errors arising from the 
stochastic method which are negligible compared to the intrinsic (gauge) noise. 
This we achieve, finding that the stochastic contribution to the total error has 
a standard deviation which is 2/3 of the standard deviation arising from the 
intrinsic variation of the signal. In the example, above we employed 24 stochastic 
sources (with no components set to zero), resulting in a cost of 24 inversions, 
per gauge configuration. Note that a similar number of inversions is needed to 
compute the (quark-connected) charged meson correlators. 
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We thus find that this variance reduction method, where applicable, is very 
powerful and effectively reduces the stochastic noise in the neutral meson corre- 
lators, making it smaller than the intrinsic noise coming from the fluctuations of 
the gauge field. 

C T-method and data-blocking 

In this appendix, we discuss the T-method and the data-blocking procedure we 
have used to estimate the statistical errors of our physical observables. 

C.l T-method 

In this section, for completeness, we just recall the basis of the T-method in- 
troduced in [21]. In the case of a primary stochastic variable with "true value" 
A (the symbol A will also be used to denote the observable itself), a suitable 
estimator of the error on the ensemble average a, i.e. its standard deviation a a , 
is given by 

1 W 

°z = x E > ( 68 ) 

n=-W 

where N is the number of measurements, 2W + 1 <C N is the number of consec- 
utive measurements used in the estimation (measurement "window") and 

JV-|n| 

r B (n) = -— V (a* - a)(a i+ H - a) . (69) 

1 1 i=l 

Here T a {n) represents the straightforward estimator of the autocorrelation func- 
tion TA(n) = ((a 1 — A)(a l+ \ n \ — A)) (the index i in a 1 labels the individual mea- 
surements, while (...) denotes the theoretical expectation value). 

The integrated autocorrelation time is conventionally defined for primary 
quantities as in eq. (125]) and estimated by (see eq. ([BHD ) 

Note that T s (0) = a^, see eq. flo^j) . is an estimate of the a priori variance of A. 

The T-method can also be applied to the analysis of secondary observables, 
F = f{A), where / denotes a non-linear function of several primary observables, 
A = {A±, A 2 , . . . }. A typical example is the case where A is given by the values of 
two-point hadron correlators at different time separations, with different smear- 
ing levels, etc., while F is a suitable estimator of the hadron mass; of course, the 

12 For a discussion of these issues, see [HS] and references therein. 
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details of the function F = f(A) depend on the specific choice of the estimator, 
e.g. on the form of the fit ansatz for the correlators and the range of time sep- 
arations employed in the fit. The main point here is that the deviation of any 
given finite-statistics estimate of F, F = f(a), from the true value f(A) can be 
approximated, in the limit of large statistics, by retaining the first term of the 
Taylor expansion of /(a) around f(A), i.e. by writing 

/(a)-/(A)^^(a a -A Q ), (71) 

a a 

where a is the index labeling the primary quantities, A a , upon which / depends. 
This remark suggests to define a new quantity, Af, which is a simple linear com- 
bination of primary quantities, and the corresponding finite-statistics estimate, 
af, via the formula 



a 



_^df(A) _ _^df(A)_ 



where a a is the ensemble average of the primary stochastic variable a a (with 
"true value" A a , as above). The variance of F — f(a) will be given by 

4 = <(/(a) -f{A)f)^({a f -A f f), (73) 

where the truncation of the Taylor series produces a relative bias 0{N~ 1 ) 
which can be neglected if the number of measurements N is sufficiently large. 
A further bias of the same order of magnitude arises from the replacement 

9 qa —> in eq. (I72*|) . which is done in practice to evaluate the first 

derivatives of / with respect to the A^s. At this point Op is estimated by the 
formula that is obtained from eq. ( 1681) by replacing r a (n) with 

N~\n\ 

= N- \n\ ^ ^ f ~ S/ ^ a / +|n ' ~ S /) ' ( 74 ) 



C.2 Binning method 

In the case where a data-blocking (also called binning) procedure is instead 
adopted to account for autocorrelations, the bin-size B plays a role similar to 
that of the window W in the T-method. The integrated autocorrelation time can 
thus be estimated, for sufficiently large values of B, by 



4(g) 
24(i) 



^(F)^^±, (75) 



where crp(B) denotes the jackknife estimate of the error on F (the mean value of 
F) that is obtained upon binning the measurements into blocks of size B. 
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C.3 Error on the error: T-method vs data-blocking 



The estimator of eq. (1681) allows to reach the optimal compromise between the 
relative statistical error on a a raising with y/W, i.e. ^statC^a)/ 1 ^ ~ ^/W/N, 
and the relative systematic error (bias) decreasing exponentially with W, i.e. 
^syst(^a) ~ l/2exp (— W/t), where r is the characteristic time of slowest expo- 
nential mode of T(n) (exponential autocorrelation time). An "optimal" value, 
W opt , to be used as upper and lower bound for the sum in eq. f[6"H]) can be ob- 
tained, e.g. by gradually increasing W and inspecting "by eye" the onset of a 
plateau for function of W, or requiring minimisation of the total error 

Stot — 4tat + 4 y st UM ■ Any va lid criterion to truncate the sum necessarily corre- 
sponds to values of W opt for which the truncation errors become comparable with 
the statistical noise level on <t s . This choice corresponds to an uncertainty on the 
error on a a decreasing like ~ 0(iV -1 / 2 ). For comparison we recall that the error 
on <7a upon use of the binning method would decrease only like ~ 0(N~ 1 ^ 3 ) |24j . 
In this case in fact the optimal choice corresponds to find a compromise between 
the relative statistical error on <t s (i.e. 5 s tat (<?a)/<?a ~ \/B/2N) which increases 
with \f~B, and the relative systematic error (bias) (i.e. 5 S yst(^a) ~ r /2-B) which 
decreases with B l . 

C.4 Further remarks 

In our error analysis carried out using the T-method, we decided to compare 
different criteria for the windowing procedure in order to test in this respect the 
robustness of our estimates. One method is given by the algorithm proposed 
in [21] which is close to optimal. A second criterion, which is slightly more 
conservative, consists in stopping the procedure as soon as f (n) becomes negative 
due to statistical fluctuations. In the 15 analysed cases (5 simulation points times 
3 quantities), no systematic trend could be detected, with the two methods giving 
in most of the cases similar results. In the cases where we cyclically vary the 
time the wall source over the lattice (see sect. I2.ip in order to restore translation 
invariance in the MC time, as required by the T-method, we average beforehand 
correlators over source cycles EE We recall that the time-slice sequences used for 
the different ensembles and the value of n = t p are specified in Table HJ 

As already mentioned, the results of the T-method have been checked against 
binning procedures. For observables that are non-linear functions of the primary 
quantities, the error estimates were obtained by combining the binning procedure 
with either bootstrap-sampling (with bin sizes B = 4,8, 16, 32 in trajectory units) 
or standard jackknife. In the latter case the optimal bin size B opt was determined 
by requiring stabilization of the estimate of the error (with B opt /T int w 10 or 
larger) . 

13 We generically find correlations between consecutive measurements taken on well separated 
time-slices (e.g. by At = 12a) to be negligible. 
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Different methods give in general comparable results. In the case of the bin- 
ning+bootstrap procedure stabilization of the error is however not always evident 
at the maximal bin size (32 in trajectory units). In particular the PCAC quark 
mass turns out to be affected by significant autocorrelations (see sect. 14.21) and 
the binning procedure seems not to be able to give reliable estimates of the error. 
In this case indeed the results lie systematically below the estimates from the 
T-method. This can be understood recalling that the T-method leads to a more 
favourable dependence upon the number of measurements in the error attributed 
to the autocorrelation time than the binning method. 

In view of these findings we have decided to use the T-method for the es- 
timates of the errors on the plaquette and am PC Ac- Also the error estimates 
for fermionic quantities (other than ampcAc) quoted in sect. E] come from this 
method. However, similar results are obtained if a binning based procedure is 
employed. 

D Details of the static potential calculation 

In this appendix we provide some details on the way we compute the static 
quark- ant iquark potential from our dynamical gauge configurations. 

D.l Improved static action 

An improvement on the signal-to-noise ratio in the measurements of the Wilson 
loop can be obtained by employing suitably smeared temporal links. This can 
be viewed as a convenient modification (or improvement) of the action for static 
quarks [66j ETJ as long as gauge invariance, cubic and parity symmetries as well 
as the local conservation of the static quark number and the static quark spin 
symmetry are preserved. Under these conditions it is still guaranteed that the 
static quark action is free from 0(a) cutoff effects [68]. The statistical improve- 
ment alluded above comes from a reduction of the noise-to-signal ratio essentially 
stemming from the fact that the modified static quark action obtained via the 
use of smeared temporal links induces a self-energy mass term with a significantly 
reduced coefficient in front of the a -1 term [HZj. For our measurements we use 
the so-called HYP-improved static quark action, which is obtained by replacing 
the temporal links U±{x,xq) in the Wilson loop by HYP-smeared links [69] 



The HYP-smearing requires the specification of three parameters a = («i, «2, a^) 
and, following ref. [67], we choose a = (1.0, 1.0,0.5) throughout our calculation. 




(76) 
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D.2 Spatial smearing 



The smoothing of the spatial links has the effect of reducing excited-state con- 
tamination in the correlation functions of the Wilson loops in the potential mea- 
surements. The operators which we measured in the simulations are constructed 
using the spatial APE smearing of ref. [70]. The smoothing procedure we use 
consists in replacing every spatial link Uj(x),j = 1, 2, 3 by itself plus a sum of its 
neighbouring spatial staples and then projecting back to the nearest element in 
the SU(3) group, i.e. we write 



Here, Psu(3)Q denotes the unique projection onto the SU(3) group element W, 
which maximises ReTrfWQ^) for any 3x3 matrix Q. The smeared and SU(3) pro- 
jected link S\Uj{x) retains all the symmetry properties of the original link Uj{x) 
under gauge transformations, charge conjugation, reflections and permutations of 
the coordinate axes. The whole set of spatially smeared links, {S\Uj(x), xeL 4 }, 
forms the spatially smeared gauge field configuration. An operator O which is 
measured on a n-times iteratively smeared gauge field configuration is called 
an operator at smearing level S n , indicated by the symbol S n O. From our 
experience a good choice is to use M = 5 different smearing levels S n , with 
n — 8, 16, 24, 32, 40, and in all cases a smearing parameter X s = 0.25. 

D.3 Static quark-ant iquark pair correlators 

The matrix of static quark-antiquark pair correlation functions, each of which 
from a technical viewpoint corresponds to a spatially smeared and temporally 
improved Wilson loop, is constructed in the following way. At fixed x we first 
form smeared string (i.e. quark-antiquark pair) operators along the three spatial 
axes, connecting x with x + ri, given by 

S n Vi(x, x + ri;x ) = 

S n Ui(x,Xo)S n Ui(x + ai,x ) . . .S n Ui(x + (r - a)i,x ), £ = 1,2,3, (78) 

and improved temporal links at fixed x, connecting x with x + 1, given by 



V^x , x + t-x) = V* YP (x , x )V* YP (x, x + a)... V* YP (x, x + {t - a)) . (79) 



S 1 U j (x) 



^SU(3) 



[Uj{x) + X s ^(UkixW^x + k)Ul(x + j) (77) 
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The smeared Wilson loop is then obtained by computing 

3 

Wim(r,t) = ^ ^ TrSiVj(x, x + ri; xq)V a {x q , x + t; x + ri) 

x,xo i=l 

S m V?(x, x + ri; x + t)V^(x , x + t;x) . (80) 

Finally we define the matrix of static quark- ant iquark pair correlators according 
to the formula 

C lm (r, t) = (Wi m (r, t)) = C ml (r, t) , (81) 

where the average is over the configurations of the ensemble. Since we have 
chosen to employ M = 5 different string operators (as discussed above) and we 
are concerned with correlators where two such operators are inserted, we end up 
with a 5 x 5 matrix of static quark-antiquark pair correlators. 
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